{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 314,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# Reference: https://jupyterbook.org/interactive/hiding.html\n",
    "# Use {hide, remove}-{input, output, cell} tags to hiding content\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "%matplotlib inline\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "from IPython.display import display\n",
    "\n",
    "sns.set()\n",
    "sns.set_context('talk')\n",
    "np.set_printoptions(threshold=20, precision=2, suppress=True)\n",
    "pd.set_option('display.max_rows', 7)\n",
    "pd.set_option('display.max_columns', 8)\n",
    "pd.set_option('precision', 2)\n",
    "# This option stops scientific notation for pandas\n",
    "# pd.set_option('display.float_format', '{:.2f}'.format)\n",
    "\n",
    "def display_df(df, rows=pd.options.display.max_rows,\n",
    "               cols=pd.options.display.max_columns):\n",
    "    with pd.option_context('display.max_rows', rows,\n",
    "                           'display.max_columns', cols):\n",
    "        display(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(sec:theory_samplingVariation)=\n",
    "# Sampling Variation\n",
    "\n",
    "Probability sampling, also known as scientific sampling, uses a chance mechanism, such as drawing indistinguishable marbles from an urn, to select a sample from the access frame. The use of randomness in the selection process enables us to calculate the variation in our sample. In particular, we can compute the chance that we draw a particular sample, which is extremely useful when considering the accuracy of our data. We demonstrate these calculations with the most basic sampling method, the Simple Random Sample (SRS), and then introduce an extension of the SRS, stratified sampling. To explain both sampling methods, we use a small population of seven elements. The population is small enough that we can list all possible samples that might result from  these sampling methods and we can compute the chance of each possible sample.\n",
    "\n",
    "__EXAMPLE: SpaceX Starship Protoypes.__\n",
    "The SpaceX Starship prototypes are called $SN1$, $SN2$, ..., where $SN$ stands for \"serial number\". In the first half of 2020, seven of these prototypes were built, and before deploying them a sample of them were pressure tested. Suppose a sample of three of the population of seven Starship prototypes for testing. (Note that while this example is artificial, much of the context is based on the actual SpaceX program and the  pressure tests made on the Starship prototypes, $SN1$, $SN2$, ..., $SN7$.)  $\\blacksquare$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Random Sample (SRS)\n",
    "\n",
    "To determine which three prototypes are to be selected for pressure testing, we use simple random sampling. We use the urn model: we write a label on each marble, place all the marbles in an urn, mix them well, and draw three without looking and without replacement between draws. All of the the possible samples we could get are listed below. \n",
    "\n",
    "$$ABC ~~ ABD ~~ ABE ~~ ABF ~~ ABG ~~ ACD ~~ ACE \\\\ ACF ~~ ACG ~~ ADE ~~ ADF ~~ ADG ~~ AEF ~~ AEG \\\\ AFG ~~ BCD ~~ BCE ~~ BCF ~~ BCG ~~ BDE ~~ BDF \\\\ BDG ~~ BEF ~~ BEG ~~BFG ~~CDE ~~ CDF ~~ CDG \\\\ CEF ~~ CEG ~~ CFG ~~ DEF ~~ DEG ~~ DFG ~~ EFG  $$\n",
    "\n",
    "We use the labels $A$, $B$, etc. rather than $SN1$, $SN2$, etc. because they are shorter and easier to distinguish from one another. There are $35$  unique samples of three from our population of seven. \n",
    "By design, each of these $35$ samples is equally likely to be chosen (the marbles are indistinguishable and well mixed) so the chance of any particular sample is the same as any other sample. This means the chance must be $1/35$,\n",
    "\n",
    "$${\\mathbb{P}}(ABC) = {\\mathbb{P}}(\\textrm{ABD}) = \\cdots = {\\mathbb{P}}(\\textrm{EFG}) = \\frac{1}{35}.$$\n",
    "\n",
    "Note that we use the special symbol ${\\mathbb{P}}$ to stand for \"probability\" or \"chance\", and we read the statement ${\\mathbb{P}}(ABC)$ as the \"chance the sample contains the marbles labeled A, B, and C.\" \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We can use the urn to answer questions about the composition of a sample. For example, to find the chance that marble $A$ is in the sample, we can add up the chance of all samples that contain $A$. There are 15 of them so the chance is:\n",
    "\n",
    "$${\\mathbb{P}}(\\textrm{A is in the sample}) = \\frac{15}{35} = \\frac{3}{7}.$$\n",
    "\n",
    "And by symmetry, we have:\n",
    "\n",
    "$${\\mathbb{P}}(\\textrm{A in sample}) = {\\mathbb{P}}(\\textrm{B in sample}) = \\cdots = {\\mathbb{P}}(\\textrm{G in sample}) = \\frac {3}{7}.$$\n",
    "\n",
    "Or, you can count them to convince yourself that each unit has an equal chance of being in the sample. \n",
    "\n",
    "We now have a more formal definition of \"representative data\" that is very useful:\n",
    "__In a SRS, every sample has the same chance of being selected.__\n",
    "\n",
    ":::{note}\n",
    "\n",
    "Many people mistakingly think that the defining property of a SRS is that every unit has an equal chance of being in the sample. However, this is not the case. A SRS of $n$ units from a population of $N$ means that every possible subset of $n$ units has the same chance of being selected.\n",
    "\n",
    ":::\n",
    "\n",
    "The stratified sample is similar to the SRS. It has the extra complication that multiple samples are drawn from disjoint subsets of the population. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stratified Sample\n",
    "\n",
    "In stratified sampling, we divide the population into non-overlapping groups, called *strata* (one group is called a *stratum* and more than one are strata), and then take a simple random sample from each.  This is like having a separate urn for each stratum and drawing marbles from each urn, independently. The strata do not have to be the same size, and we need not take the same number of units from each.\n",
    "\n",
    "In our Starship example, $SN1$, $SN2$, $SN3$ and $SN4$ were constructed differently from $SN5$, $SN6$, and $SN7$ (the later three had no flaps and no nose cone). We might consider dividing our population of seven prototypes into two strata: those with and without flaps and nose cones. Our strata would be:\n",
    "$$\\textrm{Stratum}~1: \\{A, B, C, D \\} \\\\ \\textrm{Stratum}~2: \\{ E, F, G \\}$$\n",
    "Suppose we want to take a sample of size two from the first stratum and a sample of one from the second group of prototypes.  All together, we have a sample of size three from our population, but this sampling scheme gives us fewer possible samples than the SRS.\n",
    "\n",
    "$$ABE~~ABF~~ABG~~ACE~~ACF~~ACG\\\\ADE~~ADF~~ADG~~BCE~~BCF~~BCG\\\\BDE~~BDF~~BDG~~CDE~~CDF~~CDG$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "As before, each of these samples is equally likely.\n",
    "\n",
    "$${\\mathbb{P}}\\left(ABE\\right) = {\\mathbb{P}}\\left(CDG\\right) = \\frac{1}{18}.$$\n",
    "\n",
    "However, not all triples are possible in this sampling scheme. For example,\n",
    "\n",
    "$${\\mathbb{P}}\\left(\\textrm{AEF}\\right) = 0,$$\n",
    "\n",
    "since our design dictates that we choose only one unit from the second stratum. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Again, we can compute the probability that marble $A$ is in our sample by counting up all of the occurrences of $A$ among the 18 samples:\n",
    "\n",
    "$${\\mathbb{P}}\\left(A \\textrm{ in sample}\\right) = \\frac{9}{18} = \\frac{1}{2}.$$\n",
    "\n",
    "Notice that this is the chance that $A$ is chosen from the first stratum. Since two of the four marbles are selected from the first stratum, the chance is $2/4$ or $1/2$. In this sampling scheme, not all marbles have the same chance of appearing in the sample. For example, the chance that marble $F$ is chosen is not $1/2$.\n",
    "\n",
    "$${\\mathbb{P}}\\left(F \\textrm{ in sample}\\right) = \\frac{6}{18} = {\\mathbb{P}}\\left(F \\textrm{ chosen from stratum 2} \\right)=\\frac{1}{3}.$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why would we want to use stratified sampling? \n",
    "+ It guarantees at least one element from each group.\n",
    "+ When some groups in the population are more heterogenous than others, if we take larger samples from these more variable strata, we can improve the accuracy of the sample. \n",
    "+ Stratified sampling also allows us to ensure that subgroups of the population are well-represented in the sample without using human judgement to select the individuals.\n",
    "\n",
    "Despite not all samples being possible in a stratified design, each stratum's sample is representative of its stratum.  We can use our knowledge of the sampling scheme to combine the samples from each stratum in a representative way through weighting. This topic is addressed in the exercises. \n",
    "\n",
    ":::{note}\n",
    "\n",
    "It's crucial to keep the sampling scheme in mind when we compute summary statistics, make plots, and fit models for otherwise, our analysis could be flawed. \n",
    "\n",
    ":::\n",
    "\n",
    "The simple random sample is at the core of many probability sampling schemes, not just stratified sampling. For example, most government surveys use complex sampling schemes that involve multi-stages of sampling from clusters (see the exercises) and from strata. The next section, explains how to compute the distribution of a summary statistic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sampling Distribution of a Statistic\n",
    "\n",
    "Probability sampling induces a sampling distribution on any statistic that we calculate from our sample. To explain the concept, we return to our small population of Starship prototyes, and this time give each unit a value that we want to measure. We are interested in whether or not the prototype passes a pressure test. Suppose that four of the seven will fail the test. If we take a sample and give each chosen prototype a pressure test, then we might summarize our findings with, say, the proportion of prototypes in our sample that fail the test.  Each of our possible samples (remember, there are 35) gives us a summary statistic, a sample proportion. If units $A,B,D, and F$ fail the pressure test, then for each sample, we get a sample proportion that can be $0$, $1/3$, $2/3$ and $1$. For each sample, we can calculate its corresponding proprotion. Below are a few examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|    |    |    |    |    |   \n",
    "| :---        | :----   |  :--- | :--- |  :--- |   \n",
    "| Sample      | ABC       | BCE   | BDF | CEG | \n",
    "| Proportion  | 2/3       | 1/3   | 1   | 0   | "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These samples proportions can be summarized into a probability distribution table.  There are 4 samples that give us all failed tests (sample proportion of 1). These are: $ABD$ , $ABF$, $ADF$, $BDF$, so the chance of observing a sample proportion of $1$ is $4/35$.  The probability distribution table below summaries these possible values and their chances."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Sample Proportion | No. Occurrences | Chance |\n",
    "| :---: | :---: | :---: |\n",
    "| 1 | 4 | 4/35 |\n",
    "| 2/3 | 18 | 18/35 |\n",
    "|1/3 | 12 | 12/35 |\n",
    "| 0 | 1 | 1/35 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We also can display the values in this table with a probability histogram.  And, we can find the expected outcome from the SRS, with the following reasoning: \n",
    "\n",
    "$${\\mathbb{E}}(\\textrm{sample proportion}) = 1 \\times \\frac{4}{35} + \\frac{2}{3}\\times \\frac{18}{35} + \\frac{1}{3} \\times \\frac{12}{35} + 0 \\times \\frac{1}{35}\\\\ = \\frac{20}{35} \\\\ =~ \\frac{4}{7}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we use the symbol ${\\mathbb{E}}$ to represent what we might *expect* for the sample proportion. Only one in 35 samples yiuelds a 0 for the proportion, 12 of the 35 samples have one failure (the proportion is $1/3$ for these), and so on. \n",
    "\n",
    "Further, we can find the *standard error* (we call it the *SE* for short) of the sample proportion. In some sense, this is the typical deviation of the sample proportion from the expected, $4/7$. Another name for the standard error is the *root mean square error*. This alternative name tells us how to compute it (by error we mean the difference between the sample statistic and the expected value).\n",
    "\n",
    "$${\\mathbb{SE}} = \\sqrt{(1-\\frac{4}{7})^2\\times \\frac{4}{35} + (\\frac{2}{3}-\\frac{4}{7})^2\\times \\frac{18}{35} +(\\frac{1}{3}-\\frac{4}{7})^2\\times \\frac{12}{35} +(0-\\frac{4}{7})^2\\times \\frac{1}{35} } \\\\ \\approx 0.233$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SE, also called the margin of error, indicates that even though our sample proportion has a expected value that matches $4/7$, it is likely to be around $0.23$ away from $4/7$. \n",
    "\n",
    "While these calculations are relatively straight forward, we can approximate them through a simulation study. To do this, we take samples of size $3$ from our population over and over, say 100,000 times.  For each sample, we calculate the proportion of failures. So, we have 100,000 simulated sample proportions. The probability distribution table tells us that roughly $4/35$ of these $100,000$ sample proportions should $1$, about $18/35$ of them will be $2/3$, ..., and about 2,800 of the samples will be a sample with no failures ($100,000/35=2857$). \n",
    "\n",
    "- A table of the simulated proportions should look like the probability distribution above\n",
    "- The average of the simulated proportions should be close to the expected proportion, $4/7 \\approx 0.57$ .\n",
    "- The standard deviation of the simulated proportions should be close to the standard error, about $0.23$.\n",
    "\n",
    "We show the results of such a simulation below.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating Sampling Error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__The Urn Model.__ Our urn has seven marbles, one for each protoype. Since we care only about whether the prototype fails or passes the test, we can label each marble as 'fail' or 'pass', rather than $A$ through $G$.  We create this urn as an array. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 315,
   "metadata": {},
   "outputs": [],
   "source": [
    "urn = ['fail', 'fail', 'fail', 'fail', 'pass', 'pass', 'pass']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We simulate the draw of three marbles from our urn without replacement between draws using numpy's 'random.choice' as follows. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 316,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['fail', 'fail', 'fail'], dtype='<U4')"
      ]
     },
     "execution_count": 316,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.choice(urn, size=3, replace=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can keep sampling from our urn. Below we take 10 samples from the urn: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 317,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array(['fail', 'pass', 'pass'], dtype='<U4'),\n",
       " array(['pass', 'fail', 'pass'], dtype='<U4'),\n",
       " array(['pass', 'fail', 'fail'], dtype='<U4'),\n",
       " array(['fail', 'pass', 'pass'], dtype='<U4'),\n",
       " array(['pass', 'fail', 'fail'], dtype='<U4'),\n",
       " array(['pass', 'fail', 'fail'], dtype='<U4'),\n",
       " array(['pass', 'fail', 'fail'], dtype='<U4'),\n",
       " array(['pass', 'pass', 'fail'], dtype='<U4'),\n",
       " array(['fail', 'fail', 'fail'], dtype='<U4'),\n",
       " array(['fail', 'pass', 'pass'], dtype='<U4')]"
      ]
     },
     "execution_count": 317,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[np.random.choice(urn, size = 3, replace = False) for i in range(10)] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we simply want to count the number of failures in the sample, it's easier if our urn contains 1s (for fail) and 0s (for pass) so that we can sum the results of the three draws to get the number of failures in the sample. That is, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 318,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6666666666666666"
      ]
     },
     "execution_count": 318,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "urn = [1, 1, 1, 1, 0, 0, 0]\n",
    "sum(np.random.choice(urn, size=3, replace=False))/3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For our simulation, we generate 100,000 samples, and compute the proportion of dogs in each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 319,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulations = [sum(np.random.choice(urn, size=3, replace=False)) / 3\n",
    "               for i in range(100000)] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's study these 100,000 sample proportions to find their average and standard deviation and match them against what we calculated already using the probability table."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "q1c",
     "locked": true,
     "schema_version": 2,
     "solution": false
    }
   },
   "source": [
    "### Simulation Results\n",
    "We expect the simulation results to be close to our probability calculations because we have repeated the sampling process many many times. Let's first compute the average value and the standard deviation of the 100,000 sample proportions and compare them to what the theory told us: the expected proportion is 4/7 or about 0.571, and the standard error is about 0.233."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 320,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.5708633333333333, 0.2327338708816479)"
      ]
     },
     "execution_count": 320,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(simulations), np.std(simulations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our simulation matches the theoretical values quite closely.\n",
    "We can also compare the fraction of the 100,000 values that are $0$, $1/3$, $2/3$, and $1$ and\n",
    "make a histogram. These fractions should be, approximately, $1/35$, $12/35$, $18/35$, and $4/35$, or about  $0.03$, $0.34$, $0.51$, and $0.11$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 321,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_els, counts_els = np.unique(np.array(simulations), return_counts=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 322,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.  , 0.33, 0.67, 1.  ],\n",
       "       [0.03, 0.34, 0.52, 0.11]])"
      ]
     },
     "execution_count": 322,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array((unique_els, counts_els/100000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 323,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Sample Proportion')"
      ]
     },
     "execution_count": 323,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEtCAYAAAC75j/vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABXMUlEQVR4nO3deVwU9f/A8RcsrhrgASF91VBE8UBFFDXl8FbwzDLNA8XA62teiYTH1w6vxDzCNDWvn0ead+WRdxhWGmpeaaUgKmopCoIosjC/P/jufF05XBBYgffz8eiRM/uemfdndtn3fmY+M2OmKIqCEEIIUQSYmzoBIYQQwlhStIQQQhQZUrSEEEIUGVK0hBBCFBlStIQQQhQZUrSEEEIUGVK0SoCff/6Z2rVr07x5cx4/fmzqdPKFn58fbdu2LfTtxsXFkZycnG/rS09PZ8OGDfTu3ZvGjRvTqFEjunTpwvz580lMTMy37eS3p/d/SEgItWvXLtQcjh07Ru3atQ3+q1evHs2aNaN///588803mZbZtm0btWvX5tixY7ne3rVr14yKq127NiEhIdlO54enczHV34MpWJg6AVHwdu7cyUsvvUR8fDyHDh3Cx8fH1Ck9t+HDh/Pw4cNC3WZ4eDhBQUFs376dl156KV/WGRwczO7du/H19aVbt26Ym5tz7tw5li9fzvfff8/GjRupWLFivmyrIPXp04cWLVqYZNsdOnSgQ4cOAOh0OuLi4jhw4ADBwcGcPHmSjz76SI1t2rQpoaGhODk55WobU6dOJTo6mrVr1z4zNjQ0FAcHh9w1Ihe2bt3KRx99xJkzZ9R5pvh7MBUpWsXc48eP2bdvHz169GDnzp1s3769WBQtDw+PQt/mmTNnuH//fr6t7+TJk3z33XeEhIQwePBgg9e8vb0ZO3Ysy5cvZ8KECfm2zYLi5uaGm5ubSbZdu3ZtevToYTAvMDCQ999/n40bN9K8eXM6d+4MwKuvvsqrr76a621ERERQpUoVo2KfziW//frrr6SkpBjMM8Xfg6nI4cFiLjw8nPv379O8eXM8PT358ccfuX37tqnTEsCpU6eArL9wfH19sbe357fffivkrIoHc3NzPvjgA8qXL8+XX35p6nREPpKiVcx99913mJmZ0bRpUzp06EBaWlqmY/1t27Zl8uTJbN68mXbt2tGoUSPefvttfvnllzzHTZkyhUmTJtGgQQO8vb25e/cuAJGRkfj7+6u/zAcOHMivv/6qLnvkyBFq167NmDFjDNb5n//8h9q1a3PkyBEg8zF8Pz8/hg0bxoEDB+jevTsNGjSgS5cuhIeHk5SUxNSpU2nWrBktWrRg6tSpPHr0SF1WURQ2bNhAr169cHNzo0GDBvj4+LBs2TL0dzkLCQnh888/B6Bdu3b4+fmpy1+6dImRI0fi7u6Oq6srb7/9Nj/++OMz3xtLS0sANm3aRHp6eqbXDxw4wPr16w3m/fzzzwQGBtK8eXNcXFzw8vJi6tSpBj3AkJAQunbtyokTJ+jTpw8NGzakXbt2bN++ndTUVObOnYuHhwfNmjVj7Nix3Lt3z2A/+vv7c+jQITp37kzDhg15/fXX2bt3b45tefqcVkhICD4+Ppw5c4YBAwbg6upKy5YtmT59usG+B4iKimLEiBG4u7vTvHlzpk+fzqZNm6hduzbXr19/5n7MjpWVFW3atOH333/nzp07QNbntPbu3cubb76Jm5sbTZo0YfDgwZw4cUJ9vXbt2sTGxnL8+HFq167Ntm3buH79OrVr12b16tX07duX+vXr4+/vr8ZndQ5ryZIleHl54erqysCBAw0O7+W03JPz/fz82L59e5bznz6n9ccff/Dvf/8bd3d3GjZsSO/evTlw4IBBjJ+fHwEBARw5coQ33niDBg0a0Lp1axYuXJjlZ/JFIEWrGEtKSuKHH36gUaNGvPzyy7Rq1QqtVqt+6J/0008/8fHHH9OpUyfGjBnD3bt3CQwM5Pjx43mK27VrFxcvXmTy5Mn07t0bGxsbDh48iJ+fHzdv3mTEiBGMGDGCmzdv4u/vz8GDB4GMw2I9e/bk+++/V7/4jx49yqZNm3j77bfx9vbOtr3nz59n0qRJdOzYkaCgIO7du8fYsWMZOnQosbGxjBs3jpYtW/L111+zfPlydbkFCxbw4YcfUrNmTSZOnMh7771H6dKlmTt3Ljt27AAyztnoz5tMnDiR4cOHAxlfDH369OHSpUsMGzaMcePGodPpGDp0KLt3787x/enYsSPly5dn7dq1tG/fntmzZ3PkyBF1oIdWqzWIj4iI4J133uHhw4eMHj2ayZMn07BhQ77++mtmzZplEHv79m2GDx9OkyZNeP/997GwsGDSpEkMGzaMX375hX//+9907dqVPXv2EBoaarDs5cuXGT16NE2bNiUoKAhzc3NGjx7Nd999l2N7nnb37l0CAgKoUaMGkydPpnHjxqxdu5awsDA15saNG/Tr149Tp07xzjvvEBAQwP79+5k7d26utpWdWrVqARnvU1aOHz/OuHHjsLOz4/333+fdd9/l6tWrDB48WB3sEBoaSsWKFalRowahoaE0bdpUXf6zzz7D3t6eSZMm0a1bt2zz2Lt3L6tWreLtt99m5MiRREVFMXDgQP76669ctWf48OG4u7urefXp0yfLuDNnztCnTx/OnDnD4MGDee+990hNTWXkyJGZfgj9+eefjB07lubNmzNlyhReffVVPv/8czZs2JCr3AqNIoqtLVu2KM7OzsqKFSvUeUOHDlWcnZ2V06dPq/PatGmjODs7K/v371fnxcXFKe7u7krv3r3zFFenTh0lJiZGnZeamqp4e3srrVq1UhITE9X5CQkJipeXl+Ll5aU8fvxYURRFiY+PVzw8PJQOHTood+/eVVq3bq20b99eefDggbrcgAEDlDZt2hhMOzs7K4cOHVLnrVu3TnF2djbILT09XfH29lb69OmjKIqiPH78WGncuLEybtw4g32XmJio1K9fXxk2bJg6LywsTHF2dlauXbtmsN2nc0tNTVX69euntGzZUklJSVFy8ttvvynt2rVTnJ2d1f9cXFyUYcOGGbxHiqIoAQEBSps2bTKts3fv3oqbm5s6/f777yvOzs7K2rVr1Xk//PCD4uzsnGn5t99+W/H09My0H1etWqXOe/jwodKhQwfF09NTSUtLU+Oe3P/6bT49vWbNGoNcfX19DbY3ceJEpV69esqlS5fUebdu3VIaNWqUaV8/7ZdfflGcnZ2VsLCwbGM2bdqkODs7Kzt37lQURVG2bt2qODs7K7/88ouiKIrywQcfKG5ubkp6erq6zMWLF5WOHTsqe/bsUee1adNGGTBggDp97do1xdnZWenQoYP6udVzdnZW3n//fYPpunXrKhcvXlTnXblyRXFxcVHefffdbJfLbv7T+1pRMr8fb731ltKoUSPl5s2b6rxHjx4pPXv2VBo2bKjExcWpyzk7OysHDx40iGvatKn6N/KikZ5WMab/ZazvITz576d7WzVq1KB9+/bqtI2NDT169OD06dPExcXlOs7BwcFgBNXvv//OrVu36N+/P1ZWVur8cuXKMWDAAP7++2/OnTsHQPny5fnoo4+IiYnhrbfe4tatW8yePfuZI/ZKly6Nl5eXOu3o6AhkHM7TMzMzo0qVKup5vVKlSqm9xyfdu3cPKyurHIe337t3j+PHj9OqVSsePXrE3bt3uXv3Lvfv36dDhw7cuXOHs2fP5pizq6sr33//PUuXLqVPnz5UrVqV1NRUDh8+TJ8+fQx6N0uXLmXr1q0GPbCc8nzyfa9evToAXl5eBstXrVo10zlOa2tr+vXrp06XKVOGvn378s8//6jvkbF8fX0NpuvUqaN+ThRF4eDBg3h5eRmM5rO3t6d79+652k52UlNTgYz3PSuvvPIKDx48YPr06Vy+fBnIOOy2d+9eowYsvfbaa5QqVeqZcV5eXgaHT6tVq4aXlxcRERGkpaUZ0xSj3blzh9OnT9OjRw9eeeUVdX7p0qUJCAjg0aNH/PTTT+r8smXL0rp1a4M4R0dH9ZDqi0ZGDxZT//zzD8ePH6d69eqYmZmp5wbq1KmDmZkZu3btYuLEieoXWM2aNTOto1q1aiiKQmxsLLa2trmK0/9fT799fSF5Uo0aNYCMQ0X6EWjt2rWjY8eO7Nu3j759+9K4ceNntrlChQpYWPzvI63RaLLMRaPRqOeqIKNw/fDDDxw8eJDo6GhiYmJISEgAMIh7mv7w0dq1a7MdCn3z5s1n5m1hYUHr1q3VL46oqCi++uor1q5dy/Tp0+nQoQNlypRBo9Fw7do1PvvsMy5dusTVq1f5+++/s13vk+02dl9Axg+Opw9NVqtWDYDY2FgaNmz4zDbp2djYGExrtVr1Szo+Pp74+Hi1oD5J/5l4XvHx8QDZXjYwYMAAIiIiWLduHevWraNq1aq0adOGXr16UadOnWeu/+n2ZSer9jg4OHDo0CHu3r2LnZ2dUesxRmxsLJD135r+x8GNGzfUeRUqVMDc3LD/otVqX9hzWlK0iqndu3eTlpbGlStXDHoaegkJCRw4cEAdCpzVr0X9l4v+Cy83cU/+G3L+8te/9uS6k5OT+f333wH48ccfSU5OfmZP68mC9aTsfmXrtz1hwgR27txJkyZNcHNzo0+fPjRt2pRBgwbluD19u/v372/Q+3xSVkVe7/PPP8fe3p633nrLYH6NGjWYMmUKqampbNy4kUuXLlG/fn02btzIBx98gKOjI+7u7nTs2BFXV1fWrl2b5fmmrPZHTvtCL6v3WP8F9vT7+ixPfxk+SafTAZnP3UHGr/38cOHCBczMzLK98NnKyop169bx22+/ceDAAY4cOcLatWtZv349oaGhOZ6ngtzvjycZs0/z0gvL6W9Nv80n3+Oc3qMXkRStYko/avCTTz4xOBwHcPHiRRYuXMj27dvVonX16tVM64iJiUGj0VC1alV1nrFxT9Nf4xIVFZXptejoaACDQxnz5s0jNjaW4OBg5syZw7x585gyZUpOTc6TyMhIdu7cyb///W+DEYs6nY74+Pgcr+nRt0mj0dCyZUuD1y5dusT169cpW7ZstsvrB3n06tUry2Li7OwMZBy+SUlJ4ZNPPqF58+asXLnSoCB99tlnz25oLly/fh1FUQxyunLlCvC/Hld+sLW15aWXXlLX/aSYmJjnXn9SUhIRERG4ubll2yOKjo4mMTGRRo0a0ahRI4KCgrh06RL9+/dn1apVzyxaxtL3fp4UExODtbW12gs0NzfPdMeavByiy+3fWlFTtEqsMMqVK1c4d+4czZo14/XXX6d9+/YG/w0bNgw7OzuOHj2qHl46e/aswTVBd+7c4dtvv+W1116jfPny6nxj457m4uKCnZ0dGzZsICkpSZ2flJTEV199hZ2dHfXr1wfgxIkTrF+/nt69exMQEMCbb77JunXriIyMzKc99D/6w0dP94g2bdrEw4cP1d4A/O8Xqf6XbKVKlahfvz7bt283OEyXmprKpEmTGD16tMHyT+vWrRvXrl1jyZIlmV5LSUlhx44dVK9enRo1avDo0SMePnxI9erVDQrWhQsX1JGbOW0rN+7cucOePXvU6YcPH7JhwwaqV6+er7dqMjc3p23bthw5csTgtkQJCQns3LnzudatKAozZ84kOTmZwMDAbOOmT5/Ov//9bx48eKDOq1GjBuXKlTPogZibmz/X4bIff/zR4DPy559/EhERQdu2bdUfBy+//DIXL1406CllNQJVn1d2+ej/lr799ltu3bqlzn/8+DGrVq1Cq9UW6YuRpadVDOkPFfXq1SvL10uVKsWbb77JkiVL1Gu2tFotQ4YMYdCgQZQpU4avvvqK9PR0goODDZY1Ni6rbf7nP/9h7NixvPnmm2puW7Zs4Z9//iEsLAxzc3NSUlKYPHkyNjY2BAUFARAUFMSBAweYPHky33zzDWXKlHmu/fMkNzc3rKysmDVrFjdu3KBcuXIcO3aM3bt3U7p0aYMvM/2v9eXLl+Pt7U27du2YMmUKgwYN4s0336Rv375UqFCBXbt2cfr0acaPH5/jLZiGDRvGsWPHWLBgAeHh4bRr1w4bGxtu3rzJd999x61bt1i5ciVmZmaUL18eV1dXtm3bhpWVFY6Ojvz1119s3rxZ/RJ78OBBjj8cjFWqVCkmTpzI+fPnqVSpElu3buXvv//Osrg+rzFjxhAeHk6fPn3w8/NDq9WyceNG9bozYw5n/vHHH+rnOC0tjTt37nDgwAFOnz7NwIEDszw8rjd48GCGDBlC//79ef311yldujQHDhzg6tWrzJ49W42zsbHh4sWLfPXVVzRr1izXn0GtVku/fv3w8/Pj4cOHrF69mnLlyjF27Fg1pmvXrqxcuZJ3332X1q1bc/78efbs2ZOpl6ifDgsLo3nz5lnePkv/uezVqxd9+/bF0tKSb7/9lvPnzzNlyhTKlSuXq/xfJFK0iqGdO3dibW1Nx44ds43p3bs3y5YtU0cR6m/UunjxYhITE3F3d2f8+PGZTkYbG5eVTp06sXLlShYvXsyiRYuwsLDA1dWVGTNmqNeeLFy4kOjoaObMmaP+YVWsWJEJEyYwefJkPvvsM95///287ppMXn75ZZYtW8ann37K4sWL0Wq1ODo6Mm/ePM6cOcOaNWu4c+cOL7/8Ml26dGHfvn1s27aN48eP065dO9zc3NiwYQMLFy5k1apV6HQ6HB0d+eSTT+jZs2eO2y5Tpgxr1qxhw4YN7Nmzh+XLl/PgwQNsbGxo2bIlw4YNMziZ/tlnnzFr1iy2bt3K48ePqVKlCkOHDsXJyYlRo0bxyy+/0KlTp+feJ5UqVWLSpEnMnj2b27dv4+LiwqpVqwyuT8ovDg4OrFu3jtmzZ7N06VJKly7N66+/jkajYcWKFVme73ra/v372b9/P5BRcCtVqoSjoyPz589XD39nx9PTky+++IKlS5eyePFiUlJSqFWrFvPmzaNLly5q3KhRo/jggw+YOXMmI0eOzPVhwz59+mBmZsaSJUtISUmhefPmhISEULlyZTVmzJgx6HQ6du3aRUREBK6urvzf//2f+uNNr2/fvvzyyy8sX76cs2fPZlm09J/LsLAwVq5cSXp6OnXq1GHRokXZnn8tKsyUnM7aiRKhbdu2VKlS5Zk3AzU2ThRdfn5+xMbGcujQoULZXlxcHDY2Npl6VNOmTWPDhg2cPn3aqCHlouSQc1pCCJMZM2YMXbp0MTg/8/DhQw4fPkydOnWkYIlM5PCgEMJkevTowZQpUxg6dCjt2rUjJSVFHUDw5CNFhNCToiWEMJm33nqL0qVLs2bNGubMmYO5uTn169dn9erVNGvWzNTpiReQnNMSQghRZMg5LSGEEEWGHB4sQIqikNd+rH4wVUnpB5e09kLxa7P5/Yz7NaaXy/5aseLWZmNIm3O/bE7X50nRKkCKAnFxSc8OzEL58hm3/0lIeJifKb2wSlp7ofi12a5SxoXUcf/czzamuLXZGNLm3LG1tSKna8rl8KAQQogiQ4qWEEKIIkOKlhBCiCJDzmkJIfLF7RzOZQmRX6SnJYQQosiQoiWEEKLIkKIlhMgXFdp7U6G9t6nTEMWcnNMSQuSLUmd+e2aMRmP+34dali34hLKg06Xz4EGKSbYt8ocULSFEoTEzMyM5RUd0bEKhb9uxSnm0FnJwqagzadHS6XQ0btyYlBTDXz4vvfQSp06dAiAiIoL58+dz6dIlbG1tGTBgAO+8845B/NmzZwkNDeXcuXNYWlryxhtvMGrUKINn8Vy5coVPPvmEyMhINBoNPj4+TJgwASsrKzXmzp07zJo1i4iICHQ6Ha1atWLixInY2dkV4F4QomSJjk1g0hdHC327M0d4UNuhQqFvV+Qvkxat6OhoUlJSmD17NtWrV1fnm5tn/Bo6efIkw4cPx9fXlzFjxnDixAlCQ0NRFIWAgAAAYmJi8Pf3x83NjQULFnD58mXmz59PUlISU6dOBSAhIYFBgwZhZ2fH7NmziYuLY86cOdy6dYulS5cCGQU0ICCA5ORkPvzwQ3Q6HXPnziUwMJCtW7diYSGdUiGEMDWTfhNfvHgRc3NzOnXqRNmymY9xh4WFUa9ePebMmQOAt7c3Op2OJUuW4Ofnh1arZdmyZVhbW7N48WK0Wi2tWrWiTJkyTJ8+nWHDhmFvb8/69eu5f/8+O3bsoGLFjPuj2dvbM3ToUE6fPo2rqyu7du3i4sWL7N69GycnJwDq1q1L165d2bdvH507dy68HSOEECJLJj3Ae+HCBRwcHLIsWCkpKURGRtKxY0eD+Z06deL+/fucPHkSgKNHj9KmTRu0Wq0a4+PjQ1paGhEREWpM06ZN1YIF4OnpiaWlJeHh4WpMzZo11YIFqNP6GCGEEKZl0qL1xx9/oNVqCQgIwM3NjaZNmzJ16lSSkpK4du0aqampODo6GixTrVo1IOPQ4sOHD7l582amGBsbG6ysrIiOjgYgKioqU4xGo6Fq1ao5xgA4ODioMUKI7D308+ehn7+p0xDFnMkPDyYlJfHWW28xfPhwzp07x8KFC4mOjua9994DMBgoAWBpaQlAUlISiYmJWcbo45KSMh4LkpiYaFRMzZo1s4yJiYnJU/vMzMjz0F4LCw2Q9+WLmpLWXiiGbV7+JQDZP02LHB85URgsLDSFvr+L3ftshOdp87M+IyYtWvPnz6d8+fLUrl0bgKZNm2Jra8uECRM4ejRjdFF2DwMzNzdH+e8TxrKKURRFHdCRnzFCCCFMx6RFq1mzZpnmtW7d2mBa3xN6etra2lrtPT0dA5CcnIy1tTWQ0RPLKubBgwdUqVLlmTFZ9dKMoSh5f/BbSXtwXElrLxS/NluczrhMRefqlm2MrW3e/pbyi06XVuj7u7i9z8Yolg+BjIuLY/PmzVy7ds1g/qNHjwCwtbVFo9Fw9epVg9f1046OjlhaWmJvb5/p8F1cXBxJSUnqOSpHR8dMMWlpaVy/fj3HGP32sjrXJYQwVLFDKyp2aGXqNEQxZ7KiZWZmxtSpU1m3bp3B/N27d6PRaGjZsiXu7u7s27dPPQwIsHfvXqytralfvz4AHh4eHD58mMePHxvEaDQatSfn4eHBsWPHiI+PV2MiIiJITk6mZcuWQMZowr/++ouoqCg15tKlS0RFRakxQgghTMtkhwdtbGzo378/a9euxcrKCnd3d06cOMGSJUvo378/1apVY8SIEQwePJhx48bRs2dPTp06xYoVKxg/frw6TD4wMJBdu3YxdOhQBg0axJUrV5g3bx69e/emcuXKAPTr149169bh7+/PyJEjiY+PZ86cOXh7e9O4cWMAOnfuzJIlSwgMDGT8+PEoisLcuXOpVasWvr6+ptpNQgghnmCmPNmNKWSpqamsXr2arVu3Ehsbi729Pb179yYwMFAd/LB//37CwsKIjo7G3t6e/v37Z7qNU2RkJKGhoVy4cIGKFSvy+uuvZ7qN059//snMmTM5deoUlpaWtG/fnuDgYIPzVTdv3mTGjBkcPXoUrVaLh4cHISEhVKpUKU/tS09XiIvLfJ7MGCXtOHhJay8UvzbbVSoH5PwwSFtbK85HxZn0Nk5yTqvgPe85LXPz7E9qmbRoFXdStIxX0toLxa/NUrSyVtzeZ2MUZNGSsdxCCCGKDClaQgghigy5dbkQIl/c2y/36BQFT4qWECJf5HRRsRD5RQ4PCiGEKDKkaAkh8oXV+NFYjR9t6jREMSdFSwiRL8quXU3ZtatNnYYo5qRoCSGEKDKkaAkhhCgypGgJIYQoMqRoCSGEKDKkaAkhhCgy5OJiIUS+SG3YyNQpiBJAipYQIl/EHzhi6hRECSCHB4UQQhQZUrSEEEIUGVK0hBD5wq5SOfVBkEIUFClaQgghigwpWkIIIYoMKVpCCCGKDClaQgghigwpWkIIIYqMXBWtpKQkTp06pU5HRkYyevRoxo0bR2RkZL4nJ4QQQjzJ6DtiXLp0iYEDB2Jra8t3333HtWvXGDx4MIqiUKpUKfbv38+XX35JixYtCjJfIcQLKvHTz0ydgigBjO5pLViwAIAJEyYAsHnzZnQ6HWvXruWnn36ibt26fPHFFwWSpBDixfdo4GAeDRxs6jREMWd00fr111/x9/fH29sbgEOHDlGtWjXc3NwoW7Ysr7/+OufOnSuwRIUQQgiji1ZKSgoVK1YEIDY2lkuXLuHl5WUQo9Fo8jc7IUSRUWbNKsqsWWXqNEQxZ3TRcnBw4OTJkwBs374dMzMz2rVrB4CiKHz//fdUq1atYLIUQrzwrIPGYB00xtRpiGLO6KLVt29ftm/fTrdu3fjiiy+oVasWr732Gn/++SdvvPEGkZGR+Pn55TmRd999lw4dOhjMi4iI4M0338TV1ZW2bduycuXKTMudPXsWPz8/3Nzc8PT0ZN68eaSmphrEXLlyheHDh+Pu7k7z5s354IMPSEpKMoi5c+cO48ePp3nz5jRp0oT33nuP27dv57k9Qggh8p/Rowf79u2LpaUlO3fuxM3NjZEjR6qvPXr0iGnTptGjR488JfHNN9+wf/9+HBwc1HknT55k+PDh+Pr6MmbMGE6cOEFoaCiKohAQEABATEwM/v7+uLm5sWDBAi5fvsz8+fNJSkpi6tSpACQkJDBo0CDs7OyYPXs2cXFxzJkzh1u3brF06VIAdDodAQEBJCcn8+GHH6LT6Zg7dy6BgYFs3boVCwt57JgQQrwIcvVt3L17d7p3724wz9nZmT179uQ5gb///psZM2bwyiuvGMwPCwujXr16zJkzBwBvb290Oh1LlizBz88PrVbLsmXLsLa2ZvHixWi1Wlq1akWZMmWYPn06w4YNw97envXr13P//n127NihnpOzt7dn6NChnD59GldXV3bt2sXFixfZvXs3Tk5OANStW5euXbuyb98+OnfunOf2CSGEyD+5viNGSkoKf//9Nzdu3Mjyv9yaMmUKHh4eBtd3paSkEBkZSceOHQ1iO3XqxP3799Vza0ePHqVNmzZotVo1xsfHh7S0NCIiItSYpk2bqgULwNPTE0tLS8LDw9WYmjVrqgULUKf1MUIIIUzP6J5WfHw8H330Efv37yctLS3buAsXLhi98c2bN3P+/Hl27txJaGioOv/atWukpqbi6OhoEK8f6BEdHY2rqys3b97MFGNjY4OVlRXR0dEAREVFZeodajQaqlatahDz9HogY/CJPiYvzMygfPmyeVrWwiJjJGZely9qSlp7ofi2Oaf2mJkVYiJZsLDQFPr+Lq7vc06ep83P+owYXbQ++eQT9uzZg5eXF3Xr1jXo3eRFbGwss2bNYtasWdjY2Bi8lpiYCICVlZXBfEtLSyDjdlLZxejj9AMtEhMTjYqpWbNmljExMTG5bZoQQogCYnTROnToEG+99RbTpk177o0qisKkSZNo1aoVnTp1yvJ1ALNsSq65uXmOMYqiYG7+vyOf+RWTW4oCCQkP87Ss/hdKXpcvakpae6EYtvmf+xn/z6E9traZf0AWJp0urdD3d7F7n43wPG22tbXKsbdldNHS6XQ0aNAg1wlkZf369fzxxx9899136HQ64H+FSqfTYW1tDZBpWLp+2traWu09PR0DkJycrK7Dysoqy5gHDx5QpUqVZ8Zk1UsTQghhGkZ3I5o2bcqxY8fyZaN79+7l3r17eHp64uLigouLCzt27ODq1au4uLgQGRmJRqPh6tWrBsvppx0dHbG0tMTe3j7T4bu4uDiSkpLUc1SOjo6ZYtLS0rh+/XqOMfrtZXWuSwghhGkYXbQmTZpEZGQkoaGhnDlzhuvXr+d59OBHH33Eli1bDP5r06YNr7zyClu2bMHHxwd3d3f27dun9sAgo9hZW1tTv359ADw8PDh8+DCPHz82iNFoNDRr1kyNOXbsGPHx8WpMREQEycnJtGzZEsgYTfjXX38RFRWlxly6dImoqCg1RgiRswrtvanQ3tvUaYhizujDg926dSM9PZ2VK1eyalX29xczZvRgjRo1Ms2rUKECWq1WPQQ5YsQIBg8ezLhx4+jZsyenTp1ixYoVjB8/nrJlM46XBgYGsmvXLoYOHcqgQYO4cuUK8+bNo3fv3lSuXBmAfv36sW7dOvz9/Rk5ciTx8fHMmTMHb29vGjduDEDnzp1ZsmQJgYGBjB8/HkVRmDt3LrVq1cLX19fYXSREiVbqzG+mTkGUAEYXrSFDhmQ7MKIgtGjRgoULFxIWFsbIkSOxt7cnODiYd955R41xcnJi5cqVhIaGMnr0aCpWrMjgwYMZNWqUGmNjY8OaNWuYOXMmQUFBWFpa4uPjQ3BwsBqj1WpZtWoVM2bMYMqUKWi1Wjw8PAgJCZG7YQghxAvETHny+JvIV+npCnFxmQd4GKOkjTgqae2F4tdmu0rlALitH0WYBVtbK85HxTHpi6OFlZZq5ggPajtUkNGDheB5Rw+am2ffQcp1N+LHH3/kwIED3Lhxg1KlSlG5cmVat26Np6dnrpMTQgghcsPoopWenk5QUBB79uxBURTKlStHeno6SUlJrF+/no4dO7JgwYJCPYQohBCiZDF69ODy5cvZvXs3ffv2JSIiguPHjxMZGUlERAQDBgxg7969/N///V9B5iqEEKKEM7qntW3bNtq3b68+8kPv5ZdfZvLkydy6dYstW7bg7++f3zkKIYqAh37+pk5BlABGF63Y2FgGDRqU7estWrTgyJEj+ZKUEKLoSZobZuoURAlg9OHBihUrcuXKlWxfv3LlinrrJCGEEKIgGF202rZty4YNGzh06FCm1w4ePMjGjRtp27ZtviYnhCg6LE6fwuL0KVOnIYo5ow8Pjh07lp9//pmRI0fi5OSk3pMvKiqKqKgoqlSpwtixYwsqTyHEC65ih1ZAztdpCfG8jO5pVahQgc2bNxMQEICiKBw5coTw8HDS09MZPHgwW7duzfRcLCGEECI/5eri4nLlyhEUFERQUFBB5SOEEEJkK9uidePGDWxsbChTpow6bQz9jWqFEEKI/JZt0WrXrh2hoaF069YNyBiIYczdLoy5y7sQQgiRF9kWrZEjR1K7dm2DablFkxBCCFPKtmi9++67BtNPPu4jO08+jFEIIYTIb0aPHmzXrh0HDx7M9vWdO3fi5eWVL0kJIYqee/vDubc/3NRpiGIu257W3bt3uXz5sjodGxvL2bNnKVeuXKbY9PR09u/fLz0tUSRZWpbGwsLo32/5xsJCQ3F6nJ3O1c3UKYgSINuiVbp0acaPH8/t27cBMDMzY+nSpSxdujTLeEVR6Ny5c8FkKUQBsrAw53GaQnRsQqFu17FKeV4qLU/GFiI3sv2LsbS05IsvvuDPP/9EURQmTZpE7969cXPL/GvK3NwcGxsbWrRoUaDJClFQomMTCv1pujNHeOBSw7ZQt1mQrMaPBuTGuaJg5fgzz8XFBRcXFyDjOq2OHTvi7OxcKIkJIYqWsmtXA1K0RMEy+kD+u++++8yC9fvvvz93QkIIIUR2jD6gnpqayrJly9i3bx/Jycmkp6err6WlpfHgwQOSkpLk4mIhhBAFxuie1oIFC1i4cCEJCQmULVuW2NhY/vWvf2FhYcGtW7dITU1l8uTJBZmrEEKIEs7oovX999/TrFkzDh06xJdffgnA1KlT2bt3L0uXLkWn01GqVKkCS1QIIYQwumj9/fffdOzYEXNzc+zt7bG1teXUqYwHvrVq1YqePXuyadOmAktUCCGEMLpolSlTxqAn5eDgwJ9//qlON2zYkGvXruVvdkKIIiO1YSNSGzYydRqimDO6aNWtW5cjR46o0zVq1FB7WpDRE5Mb6gpRcsUfOEL8gSPPDhTiORhdtPr378/Bgwfp168fSUlJdOnShd9//52JEyfy5Zdfsnr1aho0aFCQuQohhCjhjC5aPj4+TJs2jfj4eMqWLUvLli0ZMmQI27dvZ+7cuZQrV46JEyfmauOKorB69Wo6depEw4YN6d69O999951BTEREBG+++Saurq60bduWlStXZlrP2bNn8fPzw83NDU9PT+bNm0dqaqpBzJUrVxg+fDju7u40b96cDz74gKSkJIOYO3fuMH78eJo3b06TJk1477331NtYCSGEML1c3fjsrbfe4q233lKnx48fT9++fUlISMDJyQmtVpurjS9dupSwsDBGjRpFo0aNOHLkCEFBQWg0Gjp37szJkycZPnw4vr6+jBkzhhMnThAaGoqiKAQEBAAQExODv78/bm5uLFiwgMuXLzN//nySkpKYOnUqAAkJCQwaNAg7Oztmz55NXFwcc+bM4datW+q9FHU6HQEBASQnJ/Phhx+i0+mYO3cugYGBbN26FQsLuUecEDmxq5RxM+3b/9w3cSaiOHvub+LKlStTuXLlXC+XmprKypUr6du3LyNGjACgRYsWnDt3jnXr1tG5c2fCwsKoV68ec+bMAcDb2xudTseSJUvw8/NDq9WybNkyrK2tWbx4MVqtllatWlGmTBmmT5/OsGHDsLe3Z/369dy/f58dO3ZQsWJFAOzt7Rk6dCinT5/G1dWVXbt2cfHiRXbv3o2TkxOQcR6va9eu7Nu3T24GLIQQL4Bsi1a7du1yvTIzMzMOHDhgVKxGo2Ht2rVUqFDBYH6pUqVITk4mJSWFyMhIxo4da/B6p06dWL58OSdPnuS1117j6NGjtGnTxqCX5+Pjw0cffaQeWjx69ChNmzZVCxaAp6cnlpaWhIeH4+rqytGjR6lZs6ZasAB1Ojw8XIqWEEK8ALItWnnpPeWGubk5tWvXBjLObcXFxbFt2zZ++uknPv74Y65du0ZqaiqOjo4Gy1WrVg2A6OhoXF1duXnzZqYYGxsbrKysiI6OBiAqKoru3bsbxGg0GqpWrWoQ8/R6IGNovz4mt8zMoHz5snla1sJCA+R9+aLGlO3Vb9sUnucz8qLKqT2mHmBsYaEp9P1d0v6W4fna/KzPSLZFa+3atbneWF7t27eP0aMzHmvQunVrunfvrt7D0MrKyiDW0tISgKSkJBITE7OM0cfpB1okJiYaFVOzZs0sY2JiYvLaNCGEEPnohRhdUK9ePdatW8cff/zBZ599xtChQ9XDgtld+2Vubq4+9TWrGEVRMDf/3+DI/IrJDUWBhISHeVpW/wslr8sXNaZsryl/AT/PZ+RFY/ff/+fUHlvbzD8eC5NOl1bo+7uk/S3D87XZ1tYqx96W0UXL2HNcBw8eNHaVqldffZVXX32Vpk2bYmVlxfvvv68WpKeHpeunra2t1d7T0zEAycnJWFtbAxk9saxiHjx4QJUqVZ4Zk1UvTQghROEzumhldY4rPT2dO3fuEBMTQ/Xq1fHw8DB6w/Hx8fzwww+0aNECe3t7dX69evUAuH79OhqNhqtXrxosp592dHTE0tISe3v7TIfv4uLiSEpKUs9ROTo6ZopJS0vj+vXrdOrUSY158rZUT27P1dXV6HYJUVIlfvqZqVMQJYDRRSunc1znzp0jMDCQZs2aGb3h9PR0QkJC+Pe//62ezwI4ejTjkecNGjTA3d2dffv2MWjQIPXQ3d69e7G2tqZ+/foAeHh4cPjwYYKDg9URhHv37kWj0aj5eHh4sHLlSuLj49XRihERESQnJ9OyZUsgYzThrl27iIqKokaNGgBcunSJqKgodUi+ECJ7jwYONnUKogTI28map9SvX58BAwawaNEio5exsbGhX79+LFu2jGXLlvHzzz/z+eefM2/ePN566y1q1KjBiBEjOHnyJOPGjSM8PJwFCxawYsUKhg0bRtmyGcdMAwMDuX37NkOHDuXw4cOsWrWKWbNm0bt3b7V32K9fP7RaLf7+/uzfv5/NmzczYcIEvL29ady4MQCdO3emWrVqBAYGsmvXLnbu3MmQIUOoVasWvr6++bGbhBBCPKd8G4jx8ssvc+XKlVwtM3HiRP71r3+xZcsWFi5cyCuvvMKoUaMIDAwEMi42XrhwIWFhYYwcORJ7e3uCg4N555131HU4OTmxcuVKQkNDGT16NBUrVmTw4MGMGjVKjbGxsWHNmjXMnDmToKAgLC0t8fHxITg4WI3RarWsWrWKGTNmMGXKFLRaLR4eHoSEhMjdMIQwQpk1qwDpcYmCZaboRzw8h9u3bxMYGEhKSgrff/99fuRVLKSnK8TFZR7cYYySNuLI1KMH/7gaz6QvjhbqdmeO8MClhm2ePyMvGmNu42Rra8X5qLhC39eQsb9rO1SQ0YOF4HlHD5qbZz988LlHDz5+/Ji7d++SlpbGBx98kOsEhRBCCGM91+hByLizRPPmzenatSutW7fOr7yEEEKITPJl9KAQQghRGPJl9KAQQghRGIzuacXHxxMaGsrRo0e5ffs2WY3fMDMz4/fff8/XBIUQQgg9o4vWRx99xJ49e2jcuDHNmzdHozHdnbGFEEKUTEYXrZ9++okBAwYwZcqUgsxHCFFEyROLRWEw+pxWqVKl1NsbCSGEEKZgdNHq2bMn33zzDTqdriDzEUIIIbJl9OHBMWPGMGzYMDp16oS3tze2traZYszMzBg5cmS+JiiEKBoqtPcGIP7AERNnIoozo4vWzp07+fnnn0lPT2fDhg1ZxkjREqLkKnXmN1OnIEoAo4vW559/joODAxMnTsTR0VFGDwohhCh0Rhet27dvExISQqtWrQoyHyGEECJbRg/EqFu3LrGxsQWZixBCCJEjo4tWcHAwmzdv5quvvuKff/4hPT29IPMSQgghMjH68KD+sSPTpk1j2rRpWcbIbZyEEEIUJKOLlouLC/Xr1y/IXIQQRdhDP39TpyBKAKOL1ieffFKQeQghirikuWGmTkGUAPJoEiGEEEVGtj2tunXrEhoaSrdu3QCoU6cOZmZmOa5MzmkJUXJZnD4FgM7VzcSZiOIs26L1+uuv4+DgYDD9rKIlhCi5KnbIuIZT7vYuClK2RWvWrFkG03JOSwghhKk91zmthIQEUlJS8isXIYQQIkc5Fq3U1FQ2btzIxIkTDeZHRkbSpUsXXnvtNdzc3AgMDOTq1asFmqgQQgiRbdF6/PgxgwYN4sMPP2Tnzp3qc7SuXLlCQEAAUVFReHl54e/vT3R0NG+//TZ37twptMSFEEKUPNkWrf/7v//j1KlTTJgwgV9//RULi4zTXwsXLiQlJYUuXbqwbNkygoOD2bp1KxqNhiVLlhRa4kIIIUqebIvWnj176NSpEwEBAZQpUwbI6H0dOnQIMzMzAgIC1NgKFSrwxhtv8MMPPxR4wkIIIUqubItWTEwM7u7uBvN+++03Hj58iJ2dHXXr1jV4zcHBgX/++SdXG9c/ULJbt264ubnRvn17Zs2aRVJSkhoTERHBm2++iaurK23btmXlypWZ1nP27Fn8/Pxwc3PD09OTefPmkZqaahBz5coVhg8fjru7O82bN+eDDz4w2A7AnTt3GD9+PM2bN6dJkya899573L59O1dtEqKkurc/nHv7w02dhijmsh3ynp6enulBjz///DMALVu2zBSfmJhI2bJlc7Xx5cuXs2DBAgICAmjRogXR0dGEhYVx6dIlVqxYwcmTJxk+fDi+vr6MGTOGEydOEBoaiqIoak8vJiYGf39/3NzcWLBgAZcvX2b+/PkkJSUxdepUIGOU46BBg7Czs2P27NnExcUxZ84cbt26xdKlSwHQ6XQEBASQnJzMhx9+iE6nY+7cuQQGBrJ161b18KgQImtyUbEoDNl+Ezs4OHDhwgWDeQcOHMDMzIzWrVtnio+IiDC4GPlZFEVh+fLl9OnTh/HjxwMZxbBixYqMGzeOCxcuEBYWRr169ZgzZw4A3t7e6HQ6lixZgp+fH1qtlmXLlmFtbc3ixYvRarW0atWKMmXKMH36dIYNG4a9vT3r16/n/v377Nixg4oVKwJgb2/P0KFDOX36NK6uruzatYuLFy+ye/dunJycgIy7gnTt2pV9+/bRuXNno9smhBCiYGR7eLBLly588803HDhwgIcPH7J69Wr++usvbG1tadu2rUHst99+y9GjR2nXrp3RG37w4AHdu3ena9euBvNr1KgBwF9//UVkZCQdO3Y0eL1Tp07cv3+fkydPAnD06FHatGmDVqtVY3x8fEhLSyMiIkKNadq0qVqwADw9PbG0tCQ8PFyNqVmzplqwAHVaHyOEyJ7V+NFYjR9t6jREMZdtT8vf358ff/yRd999FzMzMxRFoVSpUsyYMUMtEPv372fdunUcP34cR0dH/P39jd6wlZUVU6ZMyTT/wIEDANSrV4/U1FQcHR0NXq9WrRoA0dHRuLq6cvPmzUwxNjY2WFlZER0dDUBUVBTdu3c3iNFoNFStWtUg5un1QEaPUx+TW2ZmUL587g6Z6llYZByazevyRY0p26vftik8z2fkRaNduxoAzfIvs40x9Z3gLCw0hb6/S9rfMjxfm5/1Gcm2aGm1WlavXs3u3bv57bffsLS0pHv37tSsWVONOXfuHCdPnqR79+6EhISoowzz6vTp0yxbtoz27duTmJgIZBS3J1laWgKQlJSUbYw+Tj/QIjEx0aiYJ9v2ZExMTMxztEoIIUR+yXF0gUajoVu3buqd3p82fPhwxowZg7n58z/h5MSJEwwfPpyqVasyffp0tXeT3U16zc3NURQl2xhFUQzyyq+Y3FAUSEh4mKdl9b9Q8rp8UWPK9pryF/DzfEZeNHb//X9O7bG1zfzjsTDpdGmFvr9L2t8yPF+bbW2tcuxtPVe1KVu2bL4UrN27dzN48GD+9a9/sXr1aipWrIi1tTVApmHp+mlra2u19/R0DEBycrK6DisrqyxjHjx4oK7DmBghhBCmZfKHQK5atYr33nuPRo0asX79eipVqgRknEvSaDSZ7mmon3Z0dMTS0hJ7e/tMh+/i4uJISkpSz1E5OjpmiklLS+P69es5xui3l9W5LiGEEIXPpEVr8+bNfPLJJ/j6+rJ8+XK1ZwRQunRp3N3d2bdvn3oYEGDv3r1YW1tTv359ADw8PDh8+DCPHz82iNFoNDRr1kyNOXbsGPHx8WpMREQEycnJ6jVnnp6e/PXXX0RFRakxly5dIioqKsvr0oQQQhQ+k10xGxcXx4wZM6hSpQr9+/fP9MRjBwcHRowYweDBgxk3bhw9e/bk1KlTrFixgvHjx6sXMgcGBrJr1y6GDh3KoEGDuHLlCvPmzaN3795UrlwZgH79+rFu3Tr8/f0ZOXIk8fHxzJkzB29vbxo3bgxA586dWbJkCYGBgYwfPx5FUZg7dy61atXC19e3cHeOEEVQasNGpk5BlADZFq0NGzbQokULqlevXiAb/vHHH3n48CGxsbH0798/0+uhoaH06NGDhQsXEhYWxsiRI7G3tyc4OJh33nlHjXNycmLlypWEhoYyevRoKlasyODBgxk1apQaY2Njw5o1a5g5cyZBQUFYWlri4+NDcHCwGqPValm1ahUzZsxgypQpaLVaPDw8CAkJkbthCGGE+ANHTJ2CKAHMlCePvT3Bzc2NyZMn06tXLwDatWvHpEmTcnUBcUmXnq4QF5d5cIcxStqII1OPHvzjajyTvjhaqNudOcIDlxq2ef6MFEW2tlacj4or9H0NGfu7tkMFGT1YCJ539KC5efbDB3O8TuvAgQM0atSIsmXLEhsby40bN7hx40aOG9QfkhNCCCHyW7ZFq1evXqxYsUK9hZGZmRkzZ85k5syZOa7w6fsVCiFKBrtK5QC4/c99E2ciirNsi9aECRNo2rQpf/zxB48fP2bRokV06NCB2rVrF2Z+QgghhCrHEQatW7dW7+i+fft2Xn/9dTmnJYQQwmSMHhZ36NAhIOOi3HPnzhEbG4tWq+WVV15Rr5kSQgghClKuxnIfPnyYjz76iL///tvgvn+VKlXigw8+yPTIEiGEECI/GV20IiMjGTVqFLa2towbNw4nJycURSEqKoqvvvqK0aNHs2bNGvViXSGEECK/GV20Fi5cSJUqVdiyZYvB7ZYg444Tb775Jl988QVffpn9s3SEEEKI52F00Tpz5gwjR47MVLAg4w7pvXr1koIlRAmW+Olnpk5BlAD5dn8iMzMzUlNT82t1Qogi5tHAwaZOQZQARt/l3dXVlS1btpCcnJzptaSkJDZv3kyDBg3yNTkhhBDiSUb3tN59910GDhxI165dGTBggHojXf1AjL///puPPvqooPIUQrzgyqxZBUiPSxQso4uWu7s7Cxcu5OOPPyY0NFR9NL2iKNjZ2TF//nxee+21AktUCPFisw4aA0jREgUrV+e02rVrR+vWrTl//jzXr18HoEqVKri4uMjjO4QQQhS4XFcajUZDw4YNadiwYUHkI4QQQmTL6IEYQgghhKlJ0RJCCFFkSNESQghRZEjREkIIUWQYXbQGDhzIzz//rE4nJSUxcOBAfv/99wJJTAhRtNz+5748tVgUuGxHD3p5eeHi4oKLiwv16tXj+PHj9O7dW309NTWV48ePk5CQUCiJCiGEENkWrYCAAC5cuMC+fftYunQpZmZmfPzxx2zatIm6devy6quvYmZmpl5kLIQQQhS0bIuWv7+/+u/Hjx/TsGFDWrdujaWlJWfOnGHLli0oisLw4cOpW7cu9evXp0GDBnTv3r0w8hZCvGAqtPcGIP7AERNnIoozoy4u1mq1QMYhw27dugFw9+5dWrZsyYABA0hLS+P8+fN88803UrSEKKFKnfnN1CmIEiDbotW7d2/q1q2Li4sLderUATA4FKj/t4eHBy1atCjgNIUQQogcilbTpk25ePEi+/fv5+7du5iZmbFgwQLCw8OpU6cOlStXlnNaQgghClW2RWvChAnqv2/dukXr1q2pVasWjx49YuPGjeoNc99//31cXV2pX78+9evXp2XLlgWftRBCiBxZWpbGwsI0l+JaWGhQFKVg1m1M0CuvvAJA586d1XNaN27coG3btnh7e/Pw4UO2bt3KggUL8nzd1oULF+jVqxcHDx5UtwcQERHB/PnzuXTpEra2tgwYMIB33nnHYNmzZ88SGhrKuXPnsLS05I033mDUqFGUKlVKjbly5QqffPIJkZGRaDQafHx8mDBhAlZWVmrMnTt3mDVrFhEREeh0Olq1asXEiROxs7PLU5uEEMJULCzMeZymEB1b+JclOVYpz0ulC+bJH0avtXLlyrz00kvqtJWVFZUrV+aNN97Azc0NyLjgOC+ioqIYNmwYOp3OYP7JkycZPnw4vr6+jBkzhhMnThAaGoqiKAQEBAAQExODv78/bm5uLFiwgMuXLzN//nySkpKYOnUqAAkJCQwaNAg7Oztmz55NXFwcc+bM4datWyxduhQAnU5HQEAAycnJfPjhh+h0OubOnUtgYCBbt26VR68IIYqc6NgEJn1xtNC3O3OEBy41bAtk3UZ/Ex86dMhguly5cpnmPdlrMYZOp+Prr79m7ty5Br0ivbCwMOrVq8ecOXMA8Pb2RqfTsWTJEvz8/NBqtSxbtgxra2sWL16MVqulVatWlClThunTpzNs2DDs7e1Zv3499+/fZ8eOHVSsWBEAe3t7hg4dyunTp3F1dWXXrl1cvHiR3bt34+TkBEDdunXp2rUr+/bto3PnzrlqmxAlzUM/f1OnIEoAk9578MSJE3z66ae88847BAUFGbyWkpJCZGQkHTt2NJjfqVMn7t+/z8mTJwE4evQobdq0UYflA/j4+JCWlkZERIQa07RpU7VgAXh6emJpaUl4eLgaU7NmTbVgAeq0PkYIkb2kuWEkzQ0zdRqimDPpMS8nJycOHDiAra0t27ZtM3jt2rVrpKam4ujoaDC/WrVqAERHR+Pq6srNmzczxdjY2GBlZUV0dDSQcfjx6evHNBoNVatWNYh5ej0ADg4OakxumZlB+fJl87SshYUGyPvyRY0p26vftik8z2ekKDL1YGMLC02h729TfbZN+bmGvH+2n/UZMWlP6+WXX8bWNuvjnomJiUDmQ46WlpZAxvmz7GL0cfpzbImJifkSI4TIntnJE5idPGHqNEQx98KOLtAPl8zuOjBzc/McYxRFwdz8fzU5v2JyQ1EgIeFhnpbV/0LJ6/JFjSnba8qezvN8Rl40di1eA8jxTu+2trk7753fdLq0Qt/fpvpsm7oHn9fPtq2tVY69rRf2eVrW1tZA5hGJ+mlra2u1Z5RVTyg5OVldh5WVVZYxDx48UNdhTIwQQgjTemGLloODAxqNhqtXrxrM1087OjpiaWmJvb09MTExBjFxcXEkJSWp56gcHR0zxaSlpXH9+vUcY/Tby+pclxBCiML3what0qVL4+7uzr59+wyurN67dy/W1tbUr18fyLj34eHDh3n8+LFBjEajoVmzZmrMsWPHiI+PV2MiIiJITk5W7+Dh6enJX3/9RVRUlBpz6dIloqKi5C4fQgjxgnhhixbAiBEjOHnyJOPGjSM8PJwFCxawYsUKhg0bRtmyGcdrAwMDuX37NkOHDuXw4cOsWrWKWbNm0bt3bypXrgxAv3790Gq1+Pv7s3//fjZv3syECRPw9vamcePGQMbdPqpVq0ZgYCC7du1i586dDBkyhFq1auHr62uyfSCEEOJ/Xuii1aJFCxYuXMjly5cZOXIk3333HcHBwQwZMkSNcXJyYuXKlSQnJzN69GhWrVrF4MGDmTx5shpjY2PDmjVrqFChAkFBQcyfPx8fHx/mz5+vxmi1WlatWkW9evWYMmUK06ZNw83NjRUrVsjdMIQQ4gXxwnwbv/HGG7zxxhuZ5nfo0IEOHTrkuKy7uzubNm3KMcbZ2ZnVq1fnGPOvf/2Lzz///Jm5CiGEMI0XpmgJIYq2e/vlzjGi4EnREkLkC52rm6lTECXAC31OSwghhHiSFC0hRL6wGj8aq/GjTZ2GKOakaAkh8kXZtaspu3a1qdMQxZwULSGEEEWGFC0hhBBFhhQtIYQQRYYULSGEEEWGFC0hhBBFhlxcLITIF6kNG5k6BVECSNESQuSL+ANHTJ2CKAHk8KAQQogiQ4qWEEKIIkOKlhAiX9hVKoddpXKmTkMUc1K0hBBCFBlStIQQQhQZUrSEEEIUGVK0hBBCFBlStIQQQhQZcnHxC0qjMcfMzIzy5csW+rZ1unQePEgp9O0KIcSzSNF6QZmZmZGcoiM6NqFQt+tYpTxaC+mAi9xL/PQzU6cgSgApWi+w6NgEJn1xtFC3OXOEB7UdKhTqNkXx8GjgYFOnIEoA+UkthBCiyJCiJYTIF2XWrKLMmlWmTkMUc3J4UAiRL6yDxgBymFAULOlpCSGEKDKkaD1l586ddOnShYYNG+Lr68uOHTtMnZIQQoj/kqL1hD179hAUFISHhweLFi2iWbNmvP/++3z//femTk0IIQRyTsvAvHnz8PX1ZdKkSQB4eXmRkJDAZ599ho+Pj4mzE0IIIT2t/7p27RpXr16lY8eOBvM7depEVFQU165dM1FmQggh9MwURVFMncSLIDw8nKFDh/LNN99Qp04ddf7vv/9Oz549+fLLL/H29s7VOvNj16anF+7bY25uVqjbe5HIvn4+ZgkZd29Rypd/Zmxh72sofvvbWEVxX5uZZb+8HB78r8TERACsrKwM5ltaWgKQlJSU63XmtOONpdGUzD80U5B9/ZwqVADAmL0o+7rwFLd9LYcH/0vfK3q60Ojnm5vLrhJCCFOTb+L/sra2BjL3qB48eGDwuhBCCNORovVfjo6OAFy9etVgfkxMjMHrQgghTEeK1n9Vq1aNqlWrZroma9++fVSvXp3KlSubKDMhhBB6MhDjCSNHjmTixImUL1+e1q1bc+jQIfbs2cP8+fNNnZoQQghkyHsmGzduZOXKldy8eZNXX32VoUOH8vrrr5s6LSGEEEjREkIIUYTIOS0hhBBFhhQtIYQQRYYULSGEEEWGFC0hhBBFhhQtIYQQRYYULRPJ7ROSHzx4wEcffYSHhwdubm4MGTKEK1euFEqu+SW3bb59+zZTpkyhTZs2uLm58cYbb7Bnz57CSTafPM+TsG/evEmTJk1YvHhxwSWYz3Lb3vT0dL744gvatWtHw4YN6datG7t27SqcZPNJbtt89+5dJk6ciKenJ82aNWPYsGFF7m9Z78KFC7i4uHDr1q0c4/L1+0sRhW737t1K7dq1lRkzZihHjhxRpk6dqjg7Oyt79uzJdpkhQ4Yor732mrJt2zZl7969Srdu3RQvLy/l/v37hZh53uW2zSkpKUr37t2VNm3aKNu2bVMiIiKU//znP4qzs7Py3XffFXL2eZOX91kvPT1d8ff3V5ydnZVFixYVQrbPLy/tnTZtmlK/fn1l5cqVyk8//aRMnjxZqV27tvLDDz8UYuZ5l9s2p6enK2+//bbSsmVLZfv27crhw4eVnj17Kl5eXkp8fHwhZ/98Ll++rHh5eSnOzs7KzZs3c4zNz+8vKVom0L59e2Xs2LEG88aMGaP4+PhkGf/rr78qzs7OSnh4uDovLi5OadSokbJ06dICzTW/5LbN+/fvV5ydnZXTp08bzA8ICFC6d+9eYHnmp9y2+Unr1q1TvL29i1TRym17Y2JilDp16iibNm0ymN+/f39l2rRpBZZnfsptm6OiohRnZ2dl+/bt6ryrV68qzs7OyrZt2woy1XyTmpqqrFu3TnFzc1OaNWv2zKKV399fcniwkOXlCclHjx7F0tISDw8PdZ6NjQ1NmzblyJEjBZ7z88pLmy0tLenTpw8NGjQwmF+jRo1MNzV+ET3Pk7CvXbvGp59+yrRp0wo6zXyTl/YeOHCAMmXKZLrjzLp165gyZUpBppsv8tLmlJQU4H/P6QMo/9+HZsbHxxdcsvnoxIkTfPrpp7zzzjsEBQU9Mz6/v7+kaBWyqKgoIPNd46tVqwZAdHR0lstUq1YNjUZjMN/BwSHL+BdNXtrcokULPv74Y4Pnm6WmphIeHk6tWrUKMNv8kZc2Q8Y5npCQEHx9fXP9pGxTykt7//jjDxwdHfnpp5/o3r079erVo2PHjuzevbvgE84HeWlznTp1aN68OYsWLeLy5cvcvXuX6dOn89JLL9G+ffuCTzofODk5ceDAAd59991M30lZye/vL7lhbiHLyxOSk5KSMsXrl8nLE5ULW349FfrTTz/lypUrLFq0KH8TLAB5bfP//d//ce3aNZYsWVKwCeazvLT37t273Lx5k0mTJjFmzBiqVq3K5s2bGTduHDY2Nrz22msFn/hzyOt7/OGHHxIYGEjnzp0B0Gq1LFq0iFdffbUAs80/L7/8cq7i8/v7S4pWIVPy8IRkJYfbQxaFJyrnpc1Px82ZM4fVq1cTEBBQJH6R5qXNUVFRLFiwgLCwsCL30NG8tDc1NZW7d++yZMkS2rRpA2T0sKOiovj8889f+KKVlzZfvnyZt99+GwcHByZNmkSZMmXYtGkTo0ePZvny5bi7uxd84oUsv7+/XvxvvGImL09ItrKyUl9/epmsfsG8aJ7nqdCPHz9m/PjxrFixgoCAAIKDgwsu0XyU2zanpaUREhKCj48PHh4e6HQ6dDodkHHIUP/vF1Ve3mNLS0s0Go3BuQ4zMzNatmzJH3/8UYDZ5o+8tHn16tUArFy5kvbt2+Pp6clnn31G3bp1mTlzZsEmbCL5/f0lRauQ5eUJyY6Ojly7di3TL5aYmJgi8UTlvD4VOikpicGDB7Nnzx4mTZpUZAoW5L7NN2/e5PTp0+zYsQMXFxf1P4CFCxeq/35R5eU9rlatWpYFOTU1NVPv5UWUlzbfuHEDJycndfAFZBTqJk2acOnSpQLM1nTy+/tLilYhy8sTkj09Pbl//z4//fSTOu/u3btERkbSsmXLAs/5eeWlzWlpaYwYMYLTp08zb948Bg0aVFjp5ovctrlSpUps2bIl038Affv2Vf/9osrLe+zl5YWiKAYXjOt0On788UeaNGlS4Dk/r7y02dHRkb/++ouEhASD+adPn6ZKlSoFmq+p5Pf3l5zTMoFnPSH57t27XL16lZo1a2JlZUXTpk1p1qwZ7733HkFBQVSoUIGFCxdibW1N3759Tdwa4+S2zRs3buT48eP06dOHf/3rX/z222/quszMzHB1dTVRS4yX2zY/Pbxfr1KlStm+9iLJbXtbtGhBq1atmD59OsnJyVSvXp2vvvqK2NhY5s6da+LWGCe3bfb39+fbb78lICCAoUOHUqZMGb755huOHz9ebJ6QXuDfX7m+skvkiw0bNigdOnRQ6tevr/j6+hpcbLh161bF2dlZ+eWXX9R58fHxSkhIiOLu7q40btxYGTJkiHL58mUTZJ53uWmzn5+f4uzsnOV/devWNVELci+37/PTitLFxYqS+/Y+fPhQ+eSTTxRPT0+lQYMGSp8+fZRjx46ZIPO8y22bL126pAwbNkxxc3NTmjRpovTt21c5evSoCTJ/fvr2PXlxcUF/f8mTi4UQQhQZck5LCCFEkSFFSwghRJEhRUsIIUSRIUVLCCFEkSFFSwghRJEhRUsIIUSRIRcXi2Lp8ePHrF69mp07d3L16lXMzc2pXr06Pj4+DBo0iNKlS5s6xSy1bduWKlWqsHbt2udeV0hICNu3bzeYZ25uTtmyZXFycqJfv3707NnzubdjCnFxcZQtW5aXXnoJ+F9bi8I9C8XzkaIlih2dTkdAQAC//fYbr7/+On369CEtLY3IyEjmzZvHoUOHWLNmDVqt1tSpFoqJEydSsWJFIOOO20lJSXz77beEhIRw79493nnnHRNnmDvh4eEEBQWxfft2tWj16dOHFi1amDgzURikaIliZ8+ePRw/fpyFCxcaPFV24MCBLF++nDlz5rBlyxb69etnwiwLT/v27alatarBvF69etG5c2cWLVrEgAEDilQBP3PmDPfv3zeY5+bmhpubm4kyEoVJzmmJYufUqVMABo+80Ovfvz+lSpUyuJdhSVSmTBnatm1LUlISf/31l6nTEcJo0tMSxY7+ybFff/11pkNfZcuW5eTJk5l6Ft9//z3r1q3jwoULpKSkUKlSJXx8fBg7dqwa6+fnx0svvcRbb71FWFgY0dHRODg4EBwcTJMmTQgNDeX7779Ho9HQoUMH9SF/kHGuqkWLFjRq1IglS5YQFxdHnTp1GDt27DMfdnjq1CnCwsLUQuvm5sbYsWNp2LDhc+0n/eM/0tLS1BxbtmxJeno63333HRUrVmTHjh3Y2NgQGRnJ559/zunTpwFo0KABo0aNomnTpur6ctNGY9f3dD516tQhPDwcgHbt2tGsWTPWrl2b5Tmt2NhYFixYwI8//siDBw9wdHRkwIAB9O7dW40JCQnht99+IzQ0lNDQUM6ePYulpSWdO3cmKChIff/Ei0N6WqLY6d69O6VKlWL27Nl07dqVBQsWcOzYMR4/fgyQqWBt3ryZMWPGYG1tTVBQEMHBwVSpUoUVK1awbNkyg9jz588zadIkOnbsSFBQEPfu3WPs2LEMHTqU2NhYxo0bR8uWLfn6669Zvny5wbI//fQTH3/8MZ06dWLMmDHcvXuXwMBAjh8/nm1bjh49ip+fH4mJiYwZM4YRI0Zw48YN+vfvT2RkZJ73UXp6OsePH0er1eLk5KTO37VrFxcvXmTy5Mn07t0bGxsbDh48iJ+fHzdv3mTEiBGMGDGCmzdv4u/vz8GDB3Pdxtys7+l8RowYQYcOHYCMc3XDhw/Psn3Xrl2jV69eHDx4kN69exMcHEz58uX5z3/+Q2hoqEHs3bt3CQgIoEaNGkyePJnGjRuzdu1awsLC8rx/RQF6vnv8CvFiOnz4sNKiRQuDu8M3atRIee+995SoqCiDWB8fH6VPnz5Kenq6Oi81NVXx9vZWunbtqs4bMGCA4uzsrBw6dEidt27dOsXZ2Vnp3bu3Oi89PV3x9vZW+vTpo85r06aN4uzsrOzfv1+dFxcXp7i7uxss26ZNG2XAgAGKoihKWlqa0q5dO+Xtt99WdDqdGvPgwQOlQ4cOSo8ePXLcB++//77i7OysnD9/XomLi1Pi4uKUf/75Rzl16pQyZswYxdnZWZk5c6bBtuvUqaPExMRk2g+tWrVSEhMT1fkJCQmKl5eX4uXlpTx+/NjoNuZ2fU/noyiKEhYWpjg7OyvXrl3L1Fa9sWPHKnXq1FHOnTunzktLS1OGDRum1K5dW/nzzz8NlluzZo3BNnx9fRVPT88c968wDelpiWKpdevWHD58mPnz59OjRw/s7OxITk5m586d9OjRw+CX/7fffsuyZcsMnpYbFxdHuXLlSE5ONlhv6dKl8fLyUqf1T15t166dOs/MzIwqVapw+/Ztg2Vr1KhB+/bt1WkbGxt69OjB6dOniYuLy9SG33//nWvXrtG+fXsSEhK4e/cud+/e5dGjR7Rp04YLFy5w69atZ+6Lnj170qJFC1q0aIGnpyd9+vRRezvjx483iHVwcMDBwcEgh1u3btG/f3+DR6OXK1eOAQMG8Pfff3Pu3Dmj25jb9T2djzHS0tL44Ycf8PT0NHjis7m5OcOHD0dRFA4dOmSwjK+vr8F0nTp1snxPhOnJOS1RbJUuXZrOnTvTuXNnIOPQ3sqVK9m5cycffPCB+sTcUqVK8euvv7Jz506ioqK4evWq+oX19NNkK1SogIXF//5sNBoNALa2tgZxGo0m0+PFa9asmSnHatWqoSgKsbGxmdahf4y7/nxLVm7evMkrr7yS436YM2cOL7/8MpDxxV2uXDmcnJyyvFbt6RyuX78OZP3o+Bo1agAZj5DXj9x7Vhtzu76n8zHGvXv3SE5OznIb+kOhsbGxBvNtbGwMprVarXquT7xYpGiJYiU5OZmlS5fi4uJiMNwdwMXFhblz53L//n2OHDnCvXv3qFixInPnzmXZsmXUq1ePRo0a0aNHD9zc3Jg2bRo3b940WMeTBetJT/bSslOqVKlM8/RfjPri96T09HQAxowZQ6NGjbJcp/6LPieNGzfONOQ9O0/n8XThzeq1J9v1rDbmdn1Z7ZdnyWkb+n369HlNc3M56FRUSNESxUrp0qVZsWIFbm5umYqWXs2aNfnxxx8pU6YMsbGxLFu2jB49emTqzdy5cydfc9P3nJ4UExODRqPJsqjoe3kvvfQSLVu2NHjtzJkzJCQkFPjoNn0OUVFRmV6Ljo4GMOjpPauNqampuVpfXtjY2PDSSy8V6DaE6cjPC1GsaDQaOnfuzPHjx/nmm28yvR4fH8/evXtp2bIlZcuWJSEhAch8WCs8PJwrV66g0+nyLbezZ88aXB92584dvv32W1577TXKly+fKb5+/frY2dmxdu1aHjx4oM5PSkpi7NixTJw4MU89kdxwcXHBzs6ODRs2kJSUZJDDV199hZ2dHfXr11fnP6uNuV1fVvS9oux6VBqNBi8vL44ePcr58+fV+Yqi8OWXX2JmZkbr1q1zsxvEC0R6WqLYCQkJ4cyZMwQHB/Ptt9/i5eWFlZUVV69eZdu2baSmpjJ16lQgo1hVrlyZJUuWkJKSwiuvvMKZM2fYvn07pUuXNigWz0ur1TJkyBAGDRpEmTJl+Oqrr0hPTyc4ODjL+FKlSvGf//yHsWPH8sYbb9CrVy9Kly7N5s2buXHjBp9++mm2hyvzy5M5vPnmm/Tq1QuALVu28M8//xAWFmZwaO1Zbczt+rKiP/+0fPlyvL29DQbB6AUFBXHs2DH8/Pzw8/PDzs6O/fv388svvzB48OAsz72JokGKlih2bGxs2LZtG6tXr+bgwYMsWrSIhw8fUqlSJTp27Mjw4cOpVKkSkPElu2zZMj755BPWrFmDoig4ODgwadIkdDodM2bM4Ny5c8/89W+MRo0a0aVLFxYvXkxiYiLu7u6MHz+eOnXqZLtMp06dWLlyJV988QWLFy/G3NycWrVq8cUXX9CmTZvnzskY+hwWL17MokWLsLCwwNXVlRkzZuDu7m4Qa0wbc7O+rHTp0oV9+/axbds2jh8/nmXRcnBwYNOmTSxYsICNGzfy6NEjnJycmDFjhlooRdFkpuR01lIIkS/y8+7tL6qS0EZhenJOSwghRJEhRUsIIUSRIUVLCCFEkSHntIQQQhQZ0tMSQghRZEjREkIIUWRI0RJCCFFkSNESQghRZEjREkIIUWRI0RJCCFFk/D8pxMoGlTHvmQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(simulations)\n",
    "plt.axvline(0.571, color='red', linestyle='dashed', linewidth=2)\n",
    "plt.title('Approximate Sampling Distribution')\n",
    "plt.ylabel('# of Simulations')\n",
    "plt.xlabel('Sample Proportion')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulations closely match the theory we developed.\n",
    "This simulation study does not *prove* the expected value of the chance outcome is 4/7 or that the chance that a sample has 2 fails is 18/35. \n",
    "But it does give us excellent approximations for these values. \n",
    "The simulation supports our earlier calculations, which is reassuring.\n",
    "More importantly, when we have a more complex setting where it might be difficult to calculate probability distributions, expected outcomes, and standard errors, a simulation study can offer valuable insights. \n",
    "\n",
    ":::{note}\n",
    "\n",
    "Simulation studies repeat a random process many many times, and summarizing patterns that result from the simulation can approximate the theoretical properties of the process. However, that is not the same as proving these theoretical properties, but often the guidance we get from the simulation is adequate for our purposes.\n",
    "\n",
    ":::\n",
    "\n",
    "Drawing marbles from an urn is such a popular framework for understanding randomness that it's probability distribution has been given a name, and most software provide the functionality to rapidly carry out the simulation.    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hypergeometric distribution\n",
    "\n",
    "The version of the urn model where we count the number of marbles of a certain type (in our case 'fail' marbles) is so common that there is a random chance process named for it: the hypergeometric. \n",
    "Instead of using `random.choice`, we can use numpy's `random.hypergeometric` to simulate drawing marbles from the urn and counting the number of `fails`. The `random.hypergeometric` method is optimzed for the 0-1 urn and allows us to ask for 100,000 simulations in one call. For completeness, we repeat our simulation study, calculate the average and standard error, and the empirical proportions.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 324,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulations_fast = np.random.hypergeometric(ngood=4, nbad=3, nsample=3, size=100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: we don't think that a pass is bad; it's just a naming convention to call the type you want to count 'good' and the other 'bad'. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 325,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5698466666666667"
      ]
     },
     "execution_count": 325,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(simulations_fast / 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 326,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.23292654741117186"
      ]
     },
     "execution_count": 326,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.std(simulations_fast / 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 327,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_els, counts_els = np.unique(np.array( simulations_fast ), return_counts=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 328,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.  , 1.  , 2.  , 3.  ],\n",
       "       [0.03, 0.35, 0.51, 0.11]])"
      ]
     },
     "execution_count": 328,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array((unique_els, counts_els/100000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perhaps the two most common probability distributions are those that arise from counting the number of 1s drawn from a 0-1 urn: drawing without replacement is the *hypergeometric* distribution and drawing with replacement is the *binomial*.  We do not delve further into the study of named probability distributions. \n",
    "\n",
    ":::{note}\n",
    "\n",
    "Whenever possible, it's a good idea to use the functionality provided in a third party package for simulating from a named distribution, rather than writing your own function, such as the random number generators offered in numpy. It's best to take advanatge of efficient and accurate code that others have  devloped.\n",
    "\n",
    ":::\n",
    "\n",
    "Our approach in this book is to develop intuition based on simulation studies to understand the results of a chance process. However, we do formalize the notion of a probability distribution, expected value, and standard deviation in Section {numref}`sec:theory_probIntro`. \n",
    "\n",
    "Now that we have simulation as a tool for understanding accuracy, we can revisit the election example from numref Chapter ch:data_scope and carry out a post-election study of what might have gone wrong with the voter polls. This simulation study imitates taking a sample of a thousand voters from a population of six million, we can examine potential sources of bias and the variation in the poll, and carry out a what-if analysis, where we examine how the predictions might have gone if even larger samples were taken.    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(sec:theory_electionpolls)=\n",
    "# Election Polls, Bias, Variance, and Big Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "load-data",
     "locked": true,
     "schema_version": 2,
     "solution": false
    }
   },
   "source": [
    "The US president is chosen by the Electoral College, not by popular vote. \n",
    "Each state is alotted a certan number of electoral college votes, as a function of their population, and\n",
    "typically, whomever wins the popular vote in the state receives all of the electoral college votes for that state. \n",
    "\n",
    "In advance of the election, polls are conducted separately in all of the states. The results help identify \"battleground\" states, and these polls are combined to predict the winner of the electoral college votes. In 2016, pollsters correctly predicted the election outcome in 46 of the 50 states. \n",
    "For those 46 states, Trump received 231 and Clinton received 232 electoral college votes.\n",
    "The remaining 4 states, Florida, Michigan, Pennsylvania, and Wisconsin, accounted for a total of 75 votes, and whichever candidate received the majority of the electoral college votes in these states would win the election. \n",
    "\n",
    "The electoral margins in these four states were narrow, e.g., in Pennsylvania Trump received 48.18% and Clinton received 47.46% of the 6,165,478 votes cast in the state. Such narrow margins can make it hard to predict the outcome given the sample sizes that the polls used.\n",
    "\n",
    "Many experts have studied the 2016 election results. According to the American Association for Public Opinion Research (AAPOR), one online, opt-in poll adjusted their polling results for education but used only three broad categories (high school or less, some college, and college graduate). They found that had they separated out those with advanced degrees from those with college degrees, then they would have reduced Clinton’s margin by 0.5 percentage points. In other words, after the fact they were able to identify an education bias where highly educated voters tended to be more willing to participate in polls. This bias matters because these voters also tended to prefer Clinton over Trump.\n",
    "\n",
    "Now that we know how people actually voted, we can carry out a simulation study that imitates the election polling under different scenarios to help develop intuition for accuracy, bias, and variance (see {cite}`grotenhuis2018`. We will simulate the polls for Pennsylvania under two scenarios: \n",
    "\n",
    "1. People surveyed didn't change their minds, didn't hide who they voted for, and were representative of those who voted on election day.\n",
    "2. People with a higher education were more likely to respond, which led to a 0.5% bias for Clinton.\n",
    "\n",
    "Our ultimate goal is to understand the chance that a poll incorrectly calls the election for Hillary Clinton even if our sample was collected with absolutely no bias, and when there is a small amount of non-response bias."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "q1",
     "locked": true,
     "schema_version": 2,
     "solution": false
    }
   },
   "source": [
    "## Simulation Study of the Sampling Error\n",
    "\n",
    "\n",
    "For the first scenario, we simulate a simple random sample (SRS) of the 6+m voters in Pennsylavania and\n",
    "calculate Trump's lead over Clinton. \n",
    "We repeat this sample collection over and over, each time calculating Trump's lead, to get \n",
    "a sense of the different values a SRS might produce.\n",
    "\n",
    "Our population consists of the votes cast for Trump, Clinton, and a third-party candidate.\n",
    "We lump all of the third party candidates together because we are only interested in the difference between\n",
    "votes cast for Trump and Clinton."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Urn Model\n",
    "\n",
    "We can cast this problem in terms of an urn model as follows:\n",
    "\n",
    "+ There are 6,165,478 marbles in the urn, one for each vote\n",
    "+ Since we care only about whether the vote is for Clinton, Trump, or some other third party candidate, we can lump all third party candidates together and label each marble in one of three ways: Trump, Clinton, and Other.\n",
    "+ The poll is a SRS of, say, 1500 marbles from the urn\n",
    "+ We tally up the counts of the three types of marbles \n",
    "\n",
    "We begin by creating an urn that represents the votes cast on election day. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 329,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2970527, 2926135,  268814])"
      ]
     },
     "execution_count": 329,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "proportions = np.array([0.4818, 0.4746, 1 - (0.4818 + 0.4746)])               \n",
    "n = 1_500\n",
    "N = 6_165_478\n",
    "votes = np.trunc(N * proportions).astype(int)\n",
    "votes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This version of the urn model is so common that there is a random chance process named for it, the multivariate hypergeometric. In Python, this is implemented in the `scipy.stats.multivariate_hypergeom.rvs` function. We can take a SRS and get our counts with the function call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 330,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([717, 725,  58])"
      ]
     },
     "execution_count": 330,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import multivariate_hypergeom\n",
    "\n",
    "multivariate_hypergeom.rvs(votes, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, each time we call `multivariate_hypergeom.rvs` we get a different sample and counts, e.g.,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 331,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([732, 700,  68])"
      ]
     },
     "execution_count": 331,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "multivariate_hypergeom.rvs(votes, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute Trumnp's lead for each simulation, i.e., $(n_T - n_C)/n$, where $n_T$ are the number of Trump votes in the sample and $n_C$ the number for Clinton. If the lead is positive, then the sample shows a win for Trump.\n",
    "\n",
    "We know the actual lead was, 0.4818 - 0.4746 =  0.0072, and to get a sense of the variation in the sampling process \n",
    "we can simulate the sampling process over and over and examine the values that we get in return. \n",
    "Below we simulate 100,000 simple random samples of 1500 voters from the state of Pennsylvania.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 332,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trump_advantage(votes, n):\n",
    "    sample_votes = multivariate_hypergeom.rvs(votes, n)\n",
    "    return (sample_votes[0] - sample_votes[1]) / n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 333,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulations = [trump_advantage(votes, n) for _ in range(100000)] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "q1c",
     "locked": true,
     "schema_version": 2,
     "solution": false
    }
   },
   "source": [
    "On average, the polling results show Trump with close to a 0.7% lead: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 334,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00723"
      ]
     },
     "execution_count": 334,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(simulations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, many times the lead in the sample was negative, meaning Clinton was the winner for that sample of voters.\n",
    "The histogram below shows the sampling distribution of Trump's advantage in Pennsylvania for a sample of 1500 voters. \n",
    "The vertical dashed line at 0 shows that more often than not, Trump is called, but there are many times when a sample \n",
    "shows Clinton in the lead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 335,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Trump Lead in the Sample')"
      ]
     },
     "execution_count": 335,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEtCAYAAAC75j/vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABQ0klEQVR4nO3deVxN+f8H8NetJFNZsjQjSuImRUpJpURabDOWEVKKUoydJMswvowh69i3ibFlmew0LSKTGYx9zUiLREaltFDd+vz+6HfPuN3Kre51u3k/H48e3M/5fM59n0+n+77nnM85Hx5jjIEQQghRAEryDoAQQgiRFCUtQgghCoOSFiGEEIVBSYsQQojCoKRFCCFEYVDSIoQQojAoaX0G/vrrLxgaGsLKygpFRUXyDkcqPD090bdv30/+vpmZmSgoKJDa+kpLSxEaGgo3NzeYm5ujW7duGDhwINatW4fc3FypvY+0le//oKAgGBoaftIYrl69CkNDQ5Gfzp07o0ePHhgzZgxOnjwp1ubYsWMwNDTE1atXq/1+qampEtUzNDREUFBQpa+loXws8vp7kAcVeQdAZO/MmTP44osvkJ2djZiYGLi6uso7pFqbOHEi3r1790nfMzY2FgEBATh+/Di++OILqawzMDAQ586dQ//+/TF48GAoKSnh/v372LVrF37//XccOnQIzZo1k8p7ydLIkSNhbW0tl/d2cnKCk5MTAEAgECAzMxPR0dEIDAzEzZs3sWTJEq6upaUlgoODYWBgUK33WLRoEZKSkrBv376P1g0ODoaurm71NqIawsLCsGTJEty9e5crk8ffg7xQ0qrnioqKEBkZiW+++QZnzpzB8ePH60XSsrW1/eTveffuXbx9+1Zq67t58yZOnz6NoKAgjBs3TmSZvb09ZsyYgV27dmHOnDlSe09ZMTMzg5mZmVze29DQEN98841Ima+vL+bOnYtDhw7BysoKAwYMAAC0bdsWbdu2rfZ7xMXFQUdHR6K65WORtr///huFhYUiZfL4e5AXOj1Yz8XGxuLt27ewsrJCr1698Mcff+D169fyDosAuHXrFoCKP3D69+8PbW1t3L59+xNHVT8oKSlh8eLFaNKkCXbu3CnvcIgUUdKq506fPg0ejwdLS0s4OTmhpKRE7Fx/3759sWDBAhw9ehSOjo7o1q0bRo0ahStXrtS43sKFCzF//nx06dIF9vb2yMrKAgBcv34d3t7e3DfzsWPH4u+//+baXrp0CYaGhpg+fbrIOr///nsYGhri0qVLAMTP4Xt6esLf3x/R0dH4+uuv0aVLFwwcOBCxsbHIy8vDokWL0KNHD1hbW2PRokV4//4915YxhtDQUHz77bcwMzNDly5d4Orqih07dkD4lLOgoCBs2rQJAODo6AhPT0+ufUJCAiZPngwLCwuYmppi1KhR+OOPPz76u1FXVwcAHDlyBKWlpWLLo6OjceDAAZGyv/76C76+vrCysoKxsTHs7OywaNEikSPAoKAgDBo0CDdu3MDIkSPRtWtXODo64vjx4yguLsaaNWtga2uLHj16YMaMGXjz5o1IP3p7eyMmJgYDBgxA165dMWTIEERERFS5LeWvaQUFBcHV1RV3796Fh4cHTE1NYWNjg2XLlon0PQAkJiZi0qRJsLCwgJWVFZYtW4YjR47A0NAQz58//2g/VkZDQwN9+vTBw4cPkZGRAaDia1oREREYPnw4zMzM0L17d4wbNw43btzglhsaGiItLQ3Xrl2DoaEhjh07hufPn8PQ0BB79uzB6NGjYWJiAm9vb65+Rdewtm3bBjs7O5iammLs2LEip/eqavdhuaenJ44fP15heflrWo8fP8Z3330HCwsLdO3aFW5uboiOjhap4+npCR8fH1y6dAnDhg1Dly5d4ODggI0bN1a4T9YFlLTqsby8PFy8eBHdunVDixYt0Lt3b6iqqnI7/Yf+/PNP/O9//4OLiwumT5+OrKws+Pr64tq1azWqd/bsWcTHx2PBggVwc3ODlpYWzp8/D09PT7x8+RKTJk3CpEmT8PLlS3h7e+P8+fMAyk6LDR06FL///jv3wX/58mUcOXIEo0aNgr29faXb++DBA8yfPx/Ozs4ICAjAmzdvMGPGDPj5+SEtLQ0zZ86EjY0NDh8+jF27dnHt1q9fjx9++AEdOnTAvHnzMGvWLDRs2BBr1qzBiRMnAJRdsxFeN5k3bx4mTpwIoOyDYeTIkUhISIC/vz9mzpwJgUAAPz8/nDt3rsrfj7OzM5o0aYJ9+/ahX79+WLlyJS5dusQN9FBVVRWpHxcXh/Hjx+Pdu3eYNm0aFixYgK5du+Lw4cP46aefROq+fv0aEydORPfu3TF37lyoqKhg/vz58Pf3x5UrV/Ddd99h0KBBCA8PR3BwsEjbp0+fYtq0abC0tERAQACUlJQwbdo0nD59usrtKS8rKws+Pj5o3749FixYAHNzc+zbtw8bNmzg6rx48QLu7u64desWxo8fDx8fH0RFRWHNmjXVeq/KdOzYEUDZ76ki165dw8yZM9GyZUvMnTsXU6ZMwbNnzzBu3DhusENwcDCaNWuG9u3bIzg4GJaWllz7n3/+Gdra2pg/fz4GDx5caRwRERHYvXs3Ro0ahcmTJyMxMRFjx47FkydPqrU9EydOhIWFBRfXyJEjK6x39+5djBw5Enfv3sW4ceMwa9YsFBcXY/LkyWJfhP755x/MmDEDVlZWWLhwIdq2bYtNmzYhNDS0WrF9MozUW7/99hvj8/nsl19+4cr8/PwYn89nd+7c4cr69OnD+Hw+i4qK4soyMzOZhYUFc3Nzq1G9Tp06sZSUFK6suLiY2dvbs969e7Pc3FyuPCcnh9nZ2TE7OztWVFTEGGMsOzub2draMicnJ5aVlcUcHBxYv379WH5+PtfOw8OD9enTR+Q1n89nMTExXNn+/fsZn88Xia20tJTZ29uzkSNHMsYYKyoqYubm5mzmzJkifZebm8tMTEyYv78/V7ZhwwbG5/NZamqqyPuWj624uJi5u7szGxsbVlhYyKpy+/Zt5ujoyPh8PvdjbGzM/P39RX5HjDHm4+PD+vTpI7ZONzc3ZmZmxr2eO3cu4/P5bN++fVzZxYsXGZ/PF2s/atQo1qtXL7F+3L17N1f27t075uTkxHr16sVKSkq4eh/2v/A9y7/eu3evSKz9+/cXeb958+axzp07s4SEBK4sPT2ddevWTayvy7ty5Qrj8/lsw4YNldY5cuQI4/P57MyZM4wxxsLCwhifz2dXrlxhjDG2ePFiZmZmxkpLS7k28fHxzNnZmYWHh3Nlffr0YR4eHtzr1NRUxufzmZOTE7ffCvH5fDZ37lyR10ZGRiw+Pp4rS05OZsbGxmzKlCmVtqusvHxfMyb++xgxYgTr1q0be/nyJVf2/v17NnToUNa1a1eWmZnJtePz+ez8+fMi9SwtLbm/kbqGjrTqMeE3Y+ERwof/L3+01b59e/Tr1497raWlhW+++QZ37txBZmZmtevp6uqKjKB6+PAh0tPTMWbMGGhoaHDljRs3hoeHB169eoX79+8DAJo0aYIlS5YgJSUFI0aMQHp6OlauXPnREXsNGzaEnZ0d91pfXx9A2ek8IR6PBx0dHe66XoMGDbijxw+9efMGGhoaVQ5vf/PmDa5du4bevXvj/fv3yMrKQlZWFt6+fQsnJydkZGTg3r17VcZsamqK33//Hdu3b8fIkSPRpk0bFBcX48KFCxg5cqTI0c327dsRFhYmcgRWVZwf/t7btWsHALCzsxNp36ZNG7FrnJqamnB3d+deq6mpYfTo0fj333+535Gk+vfvL/K6U6dO3H7CGMP58+dhZ2cnMppPW1sbX3/9dbXepzLFxcUAyn7vFfnyyy+Rn5+PZcuW4enTpwDKTrtFRERINGCpZ8+eaNCgwUfr2dnZiZw+1dPTg52dHeLi4lBSUiLJpkgsIyMDd+7cwTfffIMvv/ySK2/YsCF8fHzw/v17/Pnnn1x5o0aN4ODgIFJPX1+fO6Va19DowXrq33//xbVr19CuXTvweDzu2kCnTp3A4/Fw9uxZzJs3j/sA69Chg9g69PT0wBhDWloamjdvXq16wn+FhO8vTCQfat++PYCyU0XCEWiOjo5wdnZGZGQkRo8eDXNz849uc9OmTaGi8t8uraysXGEsysrK3LUqoCxxXbx4EefPn0dSUhJSUlKQk5MDACL1yhOePtq3b1+lQ6Ffvnz50bhVVFTg4ODAfXAkJibi4MGD2LdvH5YtWwYnJyeoqalBWVkZqamp+Pnnn5GQkIBnz57h1atXla73w+2WtC+Asi8c5U9N6unpAQDS0tLQtWvXj26TkJaWlshrVVVV7kM6Ozsb2dnZXEL9kHCfqK3s7GwAqPS2AQ8PD8TFxWH//v3Yv38/2rRpgz59+uDbb79Fp06dPrr+8ttXmYq2R1dXFzExMcjKykLLli0lWo8k0tLSAFT8tyb8cvDixQuurGnTplBSEj1+UVVVrbPXtChp1VPnzp1DSUkJkpOTRY40hHJychAdHc0NBa7o26Lww0X4gVedeh/+H6j6w1+47MN1FxQU4OHDhwCAP/74AwUFBR890vowYX2osm/ZwveeM2cOzpw5g+7du8PMzAwjR46EpaUlvLy8qnw/4XaPGTNG5OjzQxUleaFNmzZBW1sbI0aMEClv3749Fi5ciOLiYhw6dAgJCQkwMTHBoUOHsHjxYujr68PCwgLOzs4wNTXFvn37KrzeVFF/VNUXQhX9joUfYOV/rx9T/sPwQwKBAID4tTug7Nu+NDx69Ag8Hq/SG581NDSwf/9+3L59G9HR0bh06RL27duHAwcOIDg4uMrrVED1++NDkvRpTY7CqvpbE77nh7/jqn5HdRElrXpKOGpwxYoVIqfjACA+Ph4bN27E8ePHuaT17NkzsXWkpKRAWVkZbdq04cokrVee8B6XxMREsWVJSUkAIHIqY+3atUhLS0NgYCBWrVqFtWvXYuHChVVtco1cv34dZ86cwXfffScyYlEgECA7O7vKe3qE26SsrAwbGxuRZQkJCXj+/DkaNWpUaXvhII9vv/22wmTC5/MBlJ2+KSwsxIoVK2BlZYWQkBCRhPTzzz9/fEOr4fnz52CMicSUnJwM4L8jLmlo3rw5vvjiC27dH0pJSan1+vPy8hAXFwczM7NKj4iSkpKQm5uLbt26oVu3bggICEBCQgLGjBmD3bt3fzRpSUp49POhlJQUaGpqckeBSkpKYk+sqckpuur+rSkaxUqxRCLJycm4f/8+evTogSFDhqBfv34iP/7+/mjZsiUuX77MnV66d++eyD1BGRkZOHXqFHr27IkmTZpw5ZLWK8/Y2BgtW7ZEaGgo8vLyuPK8vDwcPHgQLVu2hImJCQDgxo0bOHDgANzc3ODj44Phw4dj//79uH79upR66D/C00flj4iOHDmCd+/ecUcDwH/fSIXfZFu1agUTExMcP35c5DRdcXEx5s+fj2nTpom0L2/w4MFITU3Ftm3bxJYVFhbixIkTaNeuHdq3b4/379/j3bt3aNeunUjCevToETdys6r3qo6MjAyEh4dzr9+9e4fQ0FC0a9dOqo9qUlJSQt++fXHp0iWRxxLl5OTgzJkztVo3YwzLly9HQUEBfH19K623bNkyfPfdd8jPz+fK2rdvj8aNG4scgSgpKdXqdNkff/whso/8888/iIuLQ9++fbkvBy1atEB8fLzIkVJFI1CFcVUWj/Bv6dSpU0hPT+fKi4qKsHv3bqiqqir0zch0pFUPCU8VffvttxUub9CgAYYPH45t27Zx92ypqqpiwoQJ8PLygpqaGg4ePIjS0lIEBgaKtJW0XkXv+f3332PGjBkYPnw4F9tvv/2Gf//9Fxs2bICSkhIKCwuxYMECaGlpISAgAAAQEBCA6OhoLFiwACdPnoSamlqt+udDZmZm0NDQwE8//YQXL16gcePGuHr1Ks6dO4eGDRuKfJgJv63v2rUL9vb2cHR0xMKFC+Hl5YXhw4dj9OjRaNq0Kc6ePYs7d+5g9uzZVT6Cyd/fH1evXsX69esRGxsLR0dHaGlp4eXLlzh9+jTS09MREhICHo+HJk2awNTUFMeOHYOGhgb09fXx5MkTHD16lPsQy8/Pr/KLg6QaNGiAefPm4cGDB2jVqhXCwsLw6tWrCpNrbU2fPh2xsbEYOXIkPD09oaqqikOHDnH3nUlyOvPx48fcflxSUoKMjAxER0fjzp07GDt2bIWnx4XGjRuHCRMmYMyYMRgyZAgaNmyI6OhoPHv2DCtXruTqaWlpIT4+HgcPHkSPHj2qvQ+qqqrC3d0dnp6eePfuHfbs2YPGjRtjxowZXJ1BgwYhJCQEU6ZMgYODAx48eIDw8HCxo0Th6w0bNsDKyqrCx2cJ98tvv/0Wo0ePhrq6Ok6dOoUHDx5g4cKFaNy4cbXir0soadVDZ86cgaamJpydnSut4+bmhh07dnCjCIUPat2yZQtyc3NhYWGB2bNni12MlrReRVxcXBASEoItW7Zg8+bNUFFRgampKX788Ufu3pONGzciKSkJq1at4v6wmjVrhjlz5mDBggX4+eefMXfu3Jp2jZgWLVpgx44dWL16NbZs2QJVVVXo6+tj7dq1uHv3Lvbu3YuMjAy0aNECAwcORGRkJI4dO4Zr167B0dERZmZmCA0NxcaNG7F7924IBALo6+tjxYoVGDp0aJXvraamhr179yI0NBTh4eHYtWsX8vPzoaWlBRsbG/j7+4tcTP/555/x008/ISwsDEVFRdDR0YGfnx8MDAwwdepUXLlyBS4uLrXuk1atWmH+/PlYuXIlXr9+DWNjY+zevVvk/iRp0dXVxf79+7Fy5Ups374dDRs2xJAhQ6CsrIxffvmlwutd5UVFRSEqKgpAWcJt1aoV9PX1sW7dOu70d2V69eqFrVu3Yvv27diyZQsKCwvRsWNHrF27FgMHDuTqTZ06FYsXL8by5csxefLkap82HDlyJHg8HrZt24bCwkJYWVkhKCgIrVu35upMnz4dAoEAZ8+eRVxcHExNTfHrr79yX96ERo8ejStXrmDXrl24d+9ehUlLuF9u2LABISEhKC0tRadOnbB58+ZKr78qCh6r6qod+Sz07dsXOjo6H30YqKT1iOLy9PREWloaYmJiPsn7ZWZmQktLS+yIaunSpQgNDcWdO3ckGlJOPh90TYsQIjfTp0/HwIEDRa7PvHv3DhcuXECnTp0oYRExdHqQECI333zzDRYuXAg/Pz84OjqisLCQG0Dw4ZQihAhR0iKEyM2IESPQsGFD7N27F6tWrYKSkhJMTEywZ88e9OjRQ97hkTqIrmkRQghRGHRNixBCiMKg04MyxBhDTY9jhYOp6Di4YtQ/VVN6W/bsxNLGtb9vq76ifahq8uofHq/q+/MoackQY0BmZt7HK1agSZOyx//k5LyTZkj1BvVP1Vq2KrupOfPftx+p+fmifahq8uqf5s01UNU95XR6kBBCiMKgpEUIIURhUNIihBCiMChpEVIPFRUWo6iwWN5hECJ1lLQIIYQoDEpahBBCFAYlLULqIZWePaDSkx6DROofuk+LkHpI6datKperqzeEikrNv7MKBKXIzy+scXtCaoqSFiGfIRUVJRSVMCSl5VS7rb5OE6jWIuERUhuUtAj5TCWl5WD+1svVbrd8ki0MdZtKPyBCJEBflwghhCgMSlqEEEIUBiUtQgghCoOuaRFSD5WM95F3CITIBCUtQuqhkq3byv5D026QeoaSFiEKqqp7rVRUlAH8NydSZcsJUTSUtAhRUFXda6V2/w4A4L2JaYVtO+tryTQ2QmSFkhYhCqyye61Orx0CABg860SF7UKXDZBhVITIDo0eJIQQojAoaRFCCFEYlLQIIYQoDEpahBBCFAYlLUIIIQqDkhYhhBCFQUPeCamHZoxZLe8QCJGJOnOk9ejRIxgbGyM9PV2k3MnJCYaGhmI/WVlZXJ179+7B09MTZmZm6NWrF9auXYvi4mKR9SQnJ2PixImwsLCAlZUVFi9ejLy8PJE6GRkZmD17NqysrNC9e3fMmjULr1+/lt1GEyIjT7U74Kl2B3mHQYjU1YkjrcTERPj7+0MgEIiU5+fnIzU1FbNnz0aPHj1EljVu3BgAkJKSAm9vb5iZmWH9+vV4+vQp1q1bh7y8PCxatAgAkJOTAy8vL7Rs2RIrV65EZmYmVq1ahfT0dGzfvh0AIBAI4OPjg4KCAvzwww8QCARYs2YNfH19ERYWBhWVOtFVhBDyWZPrJ7FAIMDhw4exZs0aNGjQQGz548ePwRiDo6MjDAwMKlzHjh07oKmpiS1btkBVVRW9e/eGmpoali1bBn9/f2hra+PAgQN4+/YtTpw4gWbNmgEAtLW14efnhzt37sDU1BRnz55FfHw8zp07x72XkZERBg0ahMjISAwYQE8QIIpjctRmAMBmp8lyjoQQ6ZLr6cEbN25g9erVGD9+PAICAsSWP3r0CA0bNkS7du0qXcfly5fRp08fqKqqcmWurq4oKSlBXFwcV8fS0pJLWADQq1cvqKurIzY2lqvToUMHkeQofC2sQ4iicL0XBdd7UfIOgxCpk+uRloGBAaKjo9G8eXMcO3ZMbPnjx4/RtGlTzJo1C5cvX0ZJSQkcHBwwf/58tGzZEu/evcPLly+hr68v0k5LSwsaGhpISkoCUHb68euvvxapo6ysjDZt2ojUKb8eANDV1eXqVBePV/lTtj/mY0/p/txR/8j3Se0qKsoK3/e0D1VNXv3D41W9XK5HWi1atEDz5s0rXR4fH4+MjAx07NgR27Ztw7x58/D3339j7NixeP/+PXJzcwEAGhoaYm3V1dW5gRa5ublSqUMIIUS+6vTogoULF4IxBlPTsukVLCwsYGBgAHd3d5w6dQq9e/cGAPAqSM2MMSgp/ZeTpVWnOhgDcmo4CZ/w201N29d31D/yPUIQCEoUvu9pH6qavPqneXONKo+26syQ94p07dqVS1hC3bt3h6amJuLj47kjo4qOhAoKCqCpqQmg7Eisojr5+fncOiSpQwghRL7qbNIqKChAWFgY4uPjRcoZYyguLkazZs2grq4ObW1tpKSkiNTJzMxEXl4ed41KX19frE5JSQmeP39eZR0AePbsWYXXugghhHx6dTZpNWzYECtXrsSmTZtEys+fP4/3799z923Z2triwoULKCoq4upERERAWVlZpM7Vq1eRnZ3N1YmLi0NBQQFsbGwAlI0mfPLkCRITE7k6CQkJSExM5OoQoigSWrVHQqv28g6DEKmrs0lLWVkZkyZNQlRUFJYtW4Y///wTe/bswdy5c+Ho6AgrKysAgK+vL16/fg0/Pz9cuHABu3fvxk8//QQ3Nze0bt0aAODu7g5VVVV4e3sjKioKR48exZw5c2Bvbw9zc3MAwIABA6CnpwdfX1+cPXsWZ86cwYQJE9CxY0f0799fbv1ASE3M9FiLmR5r5R0GIVJXpwdijBs3DhoaGti7dy+OHj2KJk2aYNSoUZg6dSpXx8DAACEhIQgODsa0adPQrFkzjBs3TqSOlpYW9u7di+XLlyMgIADq6upwdXVFYGAgV0dVVRW7d+/Gjz/+iIULF0JVVRW2trYICgqip2EQQkgdwWOMMXkHUV+VljJkZtZsuDyNbKoa9U9ZHzx+lo35Wy9Xu23osgFISsupUdvlk2xhqNtU4fue9qGqyXP0oJJS5cMH6+zpQUJIzZ1eOwSn1w6RdxiESB0lLUIIIQqDkhYhhBCFQUmLEEKIwqCkRQghRGFQ0iKEEKIwqpW08vLycOvWLe719evXMW3aNMycORPXr1+XenCEEELIhyS+azYhIQFjx45F8+bNcfr0aaSmpmLcuHFgjKFBgwaIiorCzp07YW1tLct4CSES2NRvkrxDIEQmJD7SWr9+PQBgzpw5AICjR49CIBBg3759+PPPP2FkZIStW7fKJEhCSPVEdHVBRFcXeYdBiNRJnLT+/vtveHt7w97eHgAQExMDPT09mJmZoVGjRhgyZAju378vs0AJIYQQiU8PFhYWolmzZgCAtLQ0JCQkwNPTU6SOsrL8pv8mhPzH5W4EAMjkaOurFupQUVGu8SSUAkEp8vMLpRwV+VxInLR0dXVx8+ZNjBgxAsePHwePx4OjoyOAsjmufv/9d+jp6cksUEKI5KZEl52ql0XSUmuogoJCAZLScqrdVl+nCVRVaNAyqTmJk9bo0aOxZMkS3L9/H4mJiejYsSN69uyJf/75B3PnzkV8fDxWrFghy1gJIXVEbR+2S0hNVStpqaur48yZMzAzM8PkyZO5Ze/fv8fSpUvxzTffyCRIQgghBKjmfFpff/01vv76a5EyPp+P8PBwqQZFCCGEVKTasxsWFhYiOzsbJSUlFS4XzhZMCCGESJvESSs7OxtLlixBVFRUpQkLAB49eiSVwAip79TVG0KlFoMSVFRotC75/EictFasWIHw8HDY2dnByMgIqqqqsoyLkHpPRUUJRSWsRqPwAKCzvpaUIyKk7pM4acXExGDEiBFYunSpLOMh5LNS01F4ABC6bEClywbPOlHDiAip2yQ+NyEQCNClSxdZxkIIIYRUSeKkZWlpiatXr8oyFkIIIaRKEp8enD9/PsaOHYvg4GC4urpCS0sLSkriOY9GDxIif+v2zwIAzPRYK+dICJEuiZPW4MGDUVpaipCQEOzevbvSejR6kBD56/BvorxDIEQmJE5aEyZMAI/Hk2UshBBCSJUkTlpTp06VZRyEEELIR1X7iRh//PEHoqOj8eLFCzRo0ACtW7eGg4MDevXqJYv4CCGEEI7ESau0tBQBAQEIDw8HYwyNGzdGaWkp8vLycODAATg7O2P9+vV0CpEQQojMSDzkfdeuXTh37hxGjx6NuLg4XLt2DdevX0dcXBw8PDwQERGBX3/9VZaxEkII+cxJfKR17Ngx9OvXD4sWLRIpb9GiBRYsWID09HT89ttv8Pb2lnaMhJBq+r2Lk7xDIEQmJE5aaWlp8PLyqnS5tbU1Ll26JJWgCCG1s9lp8scrEaKAJD492KxZMyQnJ1e6PDk5GZqamtKIiRBCCKmQxEmrb9++CA0NRUxMjNiy8+fP49ChQ+jbt69UgyOE1IzBqwQYvEqQdxiESJ3EpwdnzJiBv/76C5MnT4aBgQH09fUBAImJiUhMTISOjg5mzJghqzgJIdWw/kAAAHraO6l/JD7Satq0KY4ePQofHx8wxnDp0iXExsaitLQU48aNQ1hYGLS0aH4fQgghslOtm4sbN26MgIAABAQEyCoeQgghpFKVJq0XL15AS0sLampq3GtJ0FPeCSGEyEqlScvR0RHBwcEYPHgwgLKBGJI87YKe8k4IIURWKk1akydPhqGhochrekQTIYQQeao0aU2ZMkXktSRPeS8qKqp9RIQQQkglJB496OjoiPPnz1e6/MyZM7Czs5NKUISQ2pkxZjVmjFkt7zAIkbpKj7SysrLw9OlT7nVaWhru3buHxo0bi9UtLS1FVFQUHWkRUkc81e4g7xAIkYlKk1bDhg0xe/ZsvH79GgDA4/Gwfft2bN++vcL6jDEMGDBANlESQgghqCJpqaurY+vWrfjnn3/AGMP8+fPh5uYGMzMzsbpKSkrQ0tKCtbW1TIMlhEhmctRmAHXvwblftVCHiooymjRpVKP2AkEp8vMLpRwVUSRV3lxsbGwMY2NjAGX3aTk7O4PP53+SwAghNed6LwpA3Utaag1VUFAoQFJaTrXb6us0gaqKxJfhST0l8RMxyo8mrMjDhw/RuXPnWgVECKnfktJyMH/r5Wq3Wz7JFoa6TaUfEFEoEiet4uJi7NixA5GRkSgoKEBpaSm3rKSkBPn5+cjLy6ObiwkhhMiMxMfa69evx8aNG5GTk4NGjRohLS0NX331FVRUVJCeno7i4mIsWLBAlrESQgj5zEmctH7//Xf06NEDMTEx2LlzJwBg0aJFiIiIwPbt2yEQCNCgQQOZBUoIIYRInLRevXoFZ2dnKCkpQVtbG82bN8etW7cAAL1798bQoUNx5MiRGgfy6NEjGBsbIz09XaQ8Li4Ow4cPh6mpKfr27YuQkBCxtvfu3YOnpyfMzMzQq1cvrF27FsXFxSJ1kpOTMXHiRFhYWMDKygqLFy9GXl6eSJ2MjAzMnj0bVlZW6N69O2bNmsUN+SeEECJ/El/TUlNTEzmS0tXVxT///MO97tq1KyIiImoURGJiIvz9/SEQCETKb968iYkTJ6J///6YPn06bty4geDgYDDG4OPjAwBISUmBt7c3zMzMsH79ejx9+hTr1q1DXl4eFi1aBADIycmBl5cXWrZsiZUrVyIzMxOrVq1Ceno6d9+ZQCCAj48PCgoK8MMPP0AgEGDNmjXw9fVFWFgYVFSqNYsLIXKV0Kq9vEMgRCYk/iQ2MjLCpUuXMHLkSABA+/btuSMtoOxIrLoP1BUIBDh8+DDWrFlT4anFDRs2oHPnzli1ahUAwN7eHgKBANu2bYOnpydUVVWxY8cOaGpqYsuWLVBVVUXv3r2hpqaGZcuWwd/fH9ra2jhw4ADevn2LEydOoFmzZgAAbW1t+Pn54c6dOzA1NcXZs2cRHx+Pc+fOwcDAgNvmQYMGITIykm6cJgplpsdaeYdAiExIfHpwzJgxOH/+PNzd3ZGXl4eBAwfi4cOHmDdvHnbu3Ik9e/agS5cu1XrzGzduYPXq1Rg/frzYxJKFhYW4fv06nJ2dRcpdXFzw9u1b3Lx5EwBw+fJl9OnTB6qqqlwdV1dXlJSUIC4ujqtjaWnJJSwA6NWrF9TV1REbG8vV6dChA5ewAHCvhXUIIYTIl8RJy9XVFUuXLkV2djYaNWoEGxsbTJgwAcePH8eaNWvQuHFjzJs3r1pvbmBggOjoaEyZMgXKysoiy1JTU1FcXAx9fX2Rcj09PQBAUlIS3r17h5cvX4rV0dLSgoaGBpKSkgCUnX4sX0dZWRlt2rSpsg5QdhpUWIcQQoh8VetCzYgRIzBixAju9ezZszF69Gjk5OTAwMBA5GhHEi1atKh0WW5uLgBAQ0NDpFxdXR0AkJeXV2kdYT3hQIvc3FyJ6nToIP6QUXV1daSkpEiyOWJ4PNT4cTUqKmVJvKbt67v60D/CbZCF02uHAAAGzzohs/eQh9o8AqqidQGKvQ/Jkrz652NXmWo9uqB169Zo3bp1bVcjhjEGAJVeJ1NSUqqyDmMMSkr/HUhKqw4hhBD5qTRpOTo6VntlPB4P0dHRtQpISFNTEwDEhqULX2tqanJHT+XrAEBBQQG3Dg0NjQrr5OfnQ0dH56N1KjpKkwRjQE7Ouxq1FX67qWn7+q4+9A99w68+gaBEar/z+rAPyZK8+qd5c40qj7YqTVqyOHqqDl1dXSgrK+PZs2ci5cLX+vr6UFdXh7a2ttjpu8zMTOTl5XHXqPT19cXqlJSU4Pnz53BxceHqfDiE/8P3MzU1ldp2EUIIqblKk9a+ffs+ZRxiGjZsCAsLC0RGRsLLy4s7dRcREQFNTU2YmJgAAGxtbXHhwgUEBgZy19QiIiKgrKyMHj16cHVCQkKQnZ2Npk2bAii7abmgoAA2NjYAykYTnj17FomJiWjfvuwel4SEBCQmJmLSpEmfctMJIYRUok5frJk0aRJu3ryJmTNnIjY2FuvXr8cvv/wCf39/NGpUdujq6+uL169fw8/PDxcuXMDu3bvx008/wc3NjTtadHd3h6qqKry9vREVFYWjR49izpw5sLe3h7m5OQBgwIAB0NPTg6+vL86ePYszZ85gwoQJ6NixI/r37y+3PiCEEPIfiQdiSHqN6/z58zUOpjxra2ts3LgRGzZswOTJk6GtrY3AwECMHz+eq2NgYICQkBAEBwdj2rRpaNasGcaNG4epU6dydbS0tLB3714sX74cAQEBUFdXh6urKwIDA7k6qqqq2L17N3788UcsXLgQqqqqsLW1RVBQED0NgxBC6giJP40rusZVWlqKjIwMpKSkoF27drC1ta1xIMOGDcOwYcPEyp2cnODk5FRlWwsLi48+95DP52PPnj1V1vnqq6+wadOmj8ZKSF23qR+d0ib1k8RJq6prXPfv34evry93DYkQIl8RXV3kHQIhMiGVa1omJibw8PDA5s2bpbE6QgghpEJSG4jRokULJCcnS2t1hJBacLkbAZe7NZt1gZC6TCojDF6/fo3Q0FC539tFCCkzJXorADpNSOqfWo8eLCoqQlZWFkpKSrB48WKpBUYIIYSUV6vRg0DZ09KtrKwwaNAgODg4SCsuQgghRIxURg8SQgghn0KdfiIGIYQQ8iGJj7Sys7MRHByMy5cv4/Xr19y0IB/i8Xh4+PChVAMkhBBChCROWkuWLEF4eDjMzc1hZWUlNtMwIYQQImsSJ60///wTHh4eWLhwoSzjIYRIQX2bsZgQIYmvaTVo0ICbsoMQQgiRB4mT1tChQ3Hy5EkIBAJZxkMIIYRUSuLTg9OnT4e/vz9cXFxgb2+P5s2bi9Xh8XiYPHmyVAMkhFTfuv2zAAAzPdbKORJCpEvipHXmzBn89ddfKC0tRWhoaIV1KGmRz426ekOoqNTszhEVFdkNZurwb6LM1k2IPEmctDZt2gRdXV3MmzcP+vr6NHqQEAAqKkooKmFISsupdtvO+loyiIiQ+k3ipPX69WsEBQWhd+/esoyHEIWTlJaD+VsvV7td6LIBMoiGkPpN4vMaRkZGSEtLk2UshBBCSJUkTlqBgYE4evQoDh48iH///RelpaWyjIsQQggRI/HpQeG0I0uXLsXSpUsrrEOPcSKEECJLEictY2NjmJiYyDIWQoiU/N7FSd4hECITEietFStWyDIOQogUbXaiW09I/URTkxBCCFEYlR5pGRkZITg4GIMHDwYAdOrUCTwer8qV0TUtQuoGg1cJAICn2h3kHAkh0lVp0hoyZAh0dXVFXn8saRFC6ob1BwIA0NPeSf1TadL66aefRF7TNS1CCCHyVqtrWjk5OSgsLJRWLIQQQkiVqkxaxcXFOHToEObNmydSfv36dQwcOBA9e/aEmZkZfH198ezZM5kGSgghhFSatIqKiuDl5YUffvgBZ86c4ebRSk5Oho+PDxITE2FnZwdvb28kJSVh1KhRyMjI+GSBE0II+fxUmrR+/fVX3Lp1C3PmzMHff/8NFZWyy18bN25EYWEhBg4ciB07diAwMBBhYWFQVlbGtm3bPlnghBBCPj+VJq3w8HC4uLjAx8cHampqAMqOvmJiYsDj8eDj48PVbdq0KYYNG4aLFy/KPGBCCCGfr0pHD6akpGDYsGEiZbdv38a7d+/QqlUrGBkZiSzT1dXFv//+K5soCSHVMmPManmHQIhMVJq0SktLxSZ6/OuvvwAANjY2YvVzc3PRqFEjKYdHCKkJuqmY1FeVnh7U1dXFo0ePRMqio6PB4/Hg4OAgVj8uLk7kZmRCCCFE2io90ho4cCA2b94Me3t72Nra4vDhw3jy5AlatGiBvn37itQ9deoULl++jOnTp8s8YELIx02O2gygfj0496sW6lBRUUaTJjU/oyMQlCI/n+4tVWSVJi1vb2/88ccfmDJlCng8HhhjaNCgAX788UeoqqoCAKKiorB//35cu3YN+vr68Pb2/lRxE0Kq4HovCkD9SlpqDVVQUChAUlpOjdrr6zSBqgo9I1zRVZq0VFVVsWfPHpw7dw63b9+Guro6vv76a3To8N+58vv37+PmzZv4+uuvERQUxI0yJIQQWUhKy8H8rZdr1Hb5JFsY6jaVbkDkk6tyPi1lZWUMHjyYe9J7eRMnTsT06dOhpETfXgghhMiexJNAVoRGCxJCCPmU6BCJEEKIwqCkRQghRGHU6vQgIaRuSmjVXt4hECITlSat0NBQWFtbo127dp8wHEKINMz0WCvvEAiRiUpPDwYHB+P69evca0dHR5w/f/6TBEUIIYRUpMr7tKKjo9GtWzc0atQIaWlpePHiBV68eFHlClu3bi31IAkhhBCgiqT17bff4pdffkFsbCwAgMfjYfny5Vi+fHmVKyz/vEJCyKd3eu0QAMDgWSfkGgch0lZp0pozZw4sLS3x+PFjFBUVYfPmzXBycoKhoeGnjI8QQgjhVDl60MHBgXui+/HjxzFkyBA4Ojp+irgIIYQQMRIPeY+JiQEAlJSU4P79+0hLS4Oqqiq+/PJLmJiYyCxAQgghRKha92lduHABS5YswatXr8AYA1B2ratVq1ZYvHix2JQl0iAQCGBubo7CQtHpBL744gvcunULQNlcXuvWrUNCQgKaN28ODw8PjB8/XqT+vXv3EBwcjPv370NdXR3Dhg3D1KlT0aBBA65OcnIyVqxYgevXr0NZWRmurq6YM2cONDQ0pL5dhBBCqk/ipHX9+nVMnToVzZs3x8yZM2FgYADGGBITE3Hw4EFMmzYNe/fuhbm5uVQDTEpKQmFhIVauXClyz5jwIb03b97ExIkT0b9/f0yfPh03btxAcHAwGGPw8fEBAKSkpMDb2xtmZmZYv349nj59inXr1iEvLw+LFi0CAOTk5MDLywstW7bEypUrkZmZiVWrViE9PR3bt2+X6jYRQgipGYmT1saNG6Gjo4PffvsNmpqaIsvc3d0xfPhwbN26FTt37pRqgPHx8VBSUoKLi0uFD+jdsGEDOnfujFWrVgEA7O3tIRAIsG3bNnh6ekJVVRU7duyApqYmtmzZAlVVVfTu3RtqampYtmwZ/P39oa2tjQMHDuDt27c4ceIEmjVrBgDQ1taGn58f7ty5A1NTU6luF6k71NUbQqWG8yypqChLORpCSFUk/ku9e/cuRowYIZawAEBDQwPffvst7ty5I9XggLIh9Lq6uhUmrMLCQly/fh3Ozs4i5S4uLnj79i1u3rwJALh8+TL69OnDTV4JAK6urigpKUFcXBxXx9LSkktYANCrVy+oq6tzw/5J/aSiooSiEobHz7Kr/SM8TV7XbOo3CZv6TZJ3GIRIndSePcjj8VBcXCyt1XEeP34MVVVV+Pj44ObNm1BRUUH//v0RGBiI9PR0FBcXQ19fX6SNnp4egLJTi6ampnj58qVYHS0tLWhoaCApKQkAkJiYiK+//lqkjrKyMtq0acPVIfVXTScXDF02QAbR1F5EVxd5h0CITEictExNTfHbb7/B3d0dX3zxhciyvLw8HD16FF26dJF6gPHx8cjLy8OIESMwceJE3L9/Hxs3bkRSUhJmzZoFAGIDJdTV1bm4cnNzK6wjrJeXlwcAyM3N/Wid6uLxgCZNajbnmPC0U03b13fS7B86xff5UFFR5vYZ+hurmrz6h8erernESWvKlCkYO3YsBg0aBA8PD25QhHAgxqtXr7BkyZLaxFqhdevWoUmTJtxNzZaWlmjevDnmzJmDy5fLvhnzKtlKJSUlkVGO5THGRGZdlqQOIYrA5W4EADriIvWPxEnLwsICGzduxP/+9z8EBwdzH/CMMbRs2RLr1q1Dz549pR5gjx49xMqENzwLlT8SEr7W1NTkjp4qOloqKCjgrtFpaGhUWCc/Px86Ojo1ip0xICfnXY3aCr/d1LR9fSfN/qmP37SnRG8FQEmrPIGghNtn6G+savLqn+bNNao82qrWNS1HR0c4ODjgwYMHeP78OQBAR0cHxsbGUFGR/tRcmZmZiImJQc+ePdG2bVuu/P379wCA5s2bQ1lZGc+ePRNpJ3ytr68PdXV1aGtrIyUlRWzdeXl53LUufX19sTolJSV4/vw5XFzoD58QQuqCap/3UlZWRteuXTFgwAAMGDAApqamMklYQNnpukWLFmH//v0i5efOnYOysjJsbGxgYWGByMhIkVFcERER0NTU5J7UYWtriwsXLqCoqEikjrKyMnckZ2tri6tXryI7O5urExcXh4KCAtjY2Mhk+wghhFRPnZ65WEtLC2PGjMG+ffugoaEBCwsL3LhxA9u2bcOYMWOgp6eHSZMmYdy4cZg5cyaGDh2KW7du4ZdffsHs2bO5YfK+vr44e/Ys/Pz84OXlheTkZKxduxZubm7cVCru7u7Yv38/vL29MXnyZGRnZ2PVqlWwt7eX+g3ThBBCaqZOJy0AmDt3LrS1tREWFoYdO3ZAW1sb06ZNg6+vLwDA2toaGzduxIYNGzB58mRoa2sjMDBQ5DFOBgYGCAkJQXBwMKZNm4ZmzZph3LhxmDp1KldHS0sLe/fuxfLlyxEQEAB1dXW4uroiMDDwk28zIYSQitX5pNWgQQNMmDABEyZMqLSOk5MTnJycqlyPhYUFjhw5UmUdPp+PPXv21CRMQgghnwCN5SaEEKIwJE5aY8eOxV9//cW9zsvLw9ixY/Hw4UOZBEYIqbnBs07QrMWkXqr09KCdnR2MjY1hbGyMzp0749q1a3Bzc+OWFxcX49q1a8jJyfkkgRJCCCGVJi0fHx88evQIkZGR2L59O3g8Hv73v//hyJEjMDIyQtu2bcHj8Sp9GgUhhBAibZUmLW9vb+7/RUVF6Nq1KxwcHKCuro67d+/it99+A2MMEydOhJGREUxMTNClSxexh84SQj69dfvLnss502OtnCMhRLokGj0onNLDzs4OgwcPBgBkZWXBxsYGHh4eKCkpwYMHD3Dy5ElKWoTUAR3+TZR3CITIRKVJy83NDUZGRjA2NkanTp0AiD5QVvh/W1tbWFtbyzhMQgghpIqkZWlpifj4eERFRSErKws8Hg/r169HbGwsOnXqhNatW9M1LUIIIZ9UpUlrzpw53P/T09Ph4OCAjh074v379zh06BD3wNy5c+fC1NQUJiYmMDExoef0EUIIkRmJrml9+eWXAIABAwZw17RevHiBvn37wt7eHu/evUNYWBjWr19P920RQgiRGYkf49S6dWuRGYs1NDTQunVrDBs2DGZmZgAqnrOKEELqgq9aqNdq5mKBoBT5+YUyi49IRuKkFRMTI/K6cePGYmUVTVdPCPn0fu9S9bM4P0dqDVVQUChAUlr1H4igr9MEqir01Lu6oM4/MJcQUn2bnSbLO4Q6KSktB/O3Xq52u+WTbGGo21T6AZFqo68OhBBCFAYlLULqIYNXCTB4lSDvMAiROjo9SEg9tP5AAADQk95JvUNHWoQQQhQGJS1CCCEKg5IWIYQQhUFJixBCiMKgpEUIIURhUNIihBCiMGjIOyH10Iwxq+UdAiEyQUmL1Avq6g2hUsNnwwkfnFqfPNXuIO8QCJEJSlqkXlBRUUJRCavRw1A762vJICJCiCxQ0iL1Rk0fhhq6bIAMopGvyVGbAdCDc0n9QwMxCKmHXO9FwfVelLzDIETqKGkRQghRGJS0CCGEKAxKWoQQQhQGJS1CCCEKg5IWIYQQhUFD3gmphxJatZd3CITIBCUtQuqhmR5r5R0CITJBpwcJIYQoDEpahBBCFAadHiSkHjq9dggAYPCsE3KNo774qoU6VFSU0aRJoxq1FwhKkZ9fKOWoPk+UtAgh5CPUGqqgoFBQowcy6+s0gWoNZyAg4ihpEUKIBGr6QOblk2xhqNtU+gF9pij9E0IIURiUtAghhCgMOj1I6oTqzjwsnG1YeGG8Ps4+TAgRR0mL1Am1mXkYoNmHCflcUNIidUZNL3QD9XP24drY1G+SvEMgRCYoaRFSD0V0dZF3CITIBA3EIIQQojDoSIuQesjlbgQAOuKqC+hpGtJFSYuQemhK9FYAlLTqAnqahnRR0iKEEBmjp2lIDyWtcs6cOYOtW7ciNTUVOjo68Pf3x5AhQ+QdlkKo7r1WH6L7rAghkqCk9YHw8HAEBARg7NixsLOzQ3R0NObOnQs1NTW4urrKO7w6rzb3WtF9VoQQSVDS+sDatWvRv39/zJ8/HwBgZ2eHnJwc/Pzzz5S0JFTT0yB0nxUh4mgQhzhKWv8vNTUVz549w6xZs0TKXVxcEB4ejtTUVLRt21ZO0X06dIqPkLqDBnGI4zHGmLyDqAtiY2Ph5+eHkydPolOnTlz5w4cPMXToUOzcuRP29vbVWqeidi2Px6txW8YYSkurv91KSmXvWZO2tW2viG0/1l459y0AoESzsdTfWxHbyvO95d1WEVX1GURHWv8vNzcXAKChoSFSrq6uDgDIy8ur9jpr8+GvqHg8HpSVa77dtWlb2/aK2LbS9k2bli2T4XsrYlt5vrc8t7k+qX/HjjUkPCoqn2iE5UpK1FWEECJv9En8/zQ1NQGIH1Hl5+eLLCeEECI/lLT+n76+PgDg2bNnIuUpKSkiywkhhMgPJa3/p6enhzZt2uD3338XKY+MjES7du3QunVrOUVGCCFEiAZifGDy5MmYN28emjRpAgcHB8TExCA8PBzr1q2Td2iEEEJAQ97FHDp0CCEhIXj58iXatm0LPz8/eowTIYTUEZS0CCGEKAy6pkUIIURhUNIihBCiMChpEUIIURiUtAghhCgMSlqEEEIUBiUtOcjPz8eSJUtga2sLMzMzTJgwAcnJydVq37dvX5w8eVJsWXJyMiZOnAgLCwtYWVlh8eLFNXrYrzzVtH/u3bsHT09PmJmZoVevXli7di2Ki4tF6ixYsACGhoZiP+VvKq9rzpw5g4EDB6Jr167o378/Tpw4UWV9SfpQIBBg/fr16N27N0xNTeHu7o67d+/KbiNkTBZ9dP369Qr3F39/f9ltiIxUt38+tHLlSnh7e4uVy2UfYuSTmzBhAuvZsyc7duwYi4iIYIMHD2Z2dnbs7du3H22bm5vLPD09GZ/PZydOnBBZlp2dzezt7dnw4cNZdHQ0O3z4MLOwsGB+fn6y2hSZqEn/JCcnM3Nzc+bj48MuXrzIfvnlF2ZiYsKWLFkiUm/YsGFs+vTp7NatWyI/b968kfFW1dy5c+eYoaEh+/HHH9mlS5fYokWLGJ/PZ+Hh4ZW2kaQPf/jhB2Zqasr27dvHzp8/zzw8PJiZmRl79uzZp9gsqZJVH+3fv59169ZNbH95+vTpp9gsqalJ/wjt27eP8fl85uXlJbZMHvsQJa1P7O+//2Z8Pp/FxsZyZZmZmaxbt25s+/btVbb9448/mLOzM+vRo0eFSWvz5s2sW7duLCsriyu7ePEi4/P57Pbt29LdEBmpaf/Mnz+f9e7dmxUWFnJlBw4cYEZGRiw9PZ0xxphAIGBdu3ZlBw8elN0GyEC/fv3YjBkzRMqmT5/OXF1dK6wvSR+mpqYyIyMjkb4oLCxkDg4ObNGiRTLYCtmSRR8xxtjChQvZiBEjZBP0J1Td/mGMsfT0dDZr1izWqVMn1r17d7GkJa99iE4PfmKXL1+Guro6bG1tuTItLS1YWlri0qVLVbb19fWFiYkJdu7cWem6LS0t0axZM66sV69eUFdXR2xsrHQ2QMZq2j+XL19Gnz59oKqqypW5urqipKQEcXFxAICkpCS8f/8ehoaGstsAKRPOqO3s7CxS7uLigsTERKSmpoq1kaQPr1y5gpKSEri4uHB1VFVV4eDg8NH9sK6RVR8BwKNHjxRqf6lITfoHANatW4eHDx9i9+7dMDIyElsur32IktYnlpiYCD09PSgri07Pp6uri6SkpCrbnjp1CmvWrIGWllal6y7/NHplZWW0adPmo+uuK2rSP+/evcPLly/Ftl1LSwsaGhpcu/j4eADAiRMn0KtXL5iYmNT56ziJiYkAxGcZ0NPTA4AK+0SSPkxMTESTJk3E9iU9PT28ePEC79+/l9o2yJqs+qi0tBRPnjxBeno6hg4dChMTEzg4OCAkJEShZiWvSf8AZV+Sz549i549e1a6XnnsQ/TAXCkSCAQ4e/ZspctbtGiBvLw8sdmRgbIZkj82YILP51e5PDc3t8br/hRk1T+VzTpdvp0waeXm5mL16tV4+/Yttm/fjrFjx+LIkSMf7V95qMmM2pL0YVV1gLJBCmpqarUL/hORVR8Jj8yTkpIwa9YsNGvWDOfPn0dwcDDy8vIwbdo0aW+KTNR0VvYOHTpUuV557UOUtKSosLAQgYGBlS7v0aMHGjRoUOlyacyOXH7mZaBs9uW6MPOyrPqHVTLrtHCZsN2IESNgaWmJ3r17c8t79uwJZ2dnbN++HWvWrJFoOz6lyraNVTGjdlVHAcL6ldWpqi/rKln1kba2Nnbu3AkjIyO0bNkSAGBtbY33799j586dGD9+fIUf2nVNTfqnOuuV9P2khZKWFKmrq+Px48dV1pk2bRqeP38uVp6fn1/rPwANDY0KvzXl5+dDR0enVuuWBln1j7C8om0vKCjgZp3W09PjTokINW7cGObm5h+NS15qMqO2hobGR/tQQ0ODW0dF61WED2MhWfaRvb29WB0HBwccPXoUSUlJ6NKlS63jlzVZzcour31I/l+/PzP6+vpITU0V+5aSkpJS69mR9fX1uZmWhUpKSvD8+XOFmXm5Jv2jrq4ObW1tsW3PzMxEXl4e1y4yMrLCASmFhYUig1fqkprMqC1JH7Zv3x7Z2dnIyckRq9OmTRuRAS11naz66PHjxzh48KDYvX7CazV1dZ8pT1azsstrH6Kk9Yn16tULb9++xZ9//smVZWVl4fr167CxsanVum1tbXH16lVkZ2dzZXFxcSgoKKj1uj+VmvaPra0tLly4gKKiIq4sIiICysrK6NGjBwDg2LFjWLhwocgF4levXuHmzZtcnbqmJjNqS9KHwn8jIiK4OkVFRYiNjVWYfUVIVn2UkpKCJUuWiI2EO3fuHNq0aVMnzl5IQlazssttH5LZYHpSKQ8PD9ajRw925MgRFhkZyd3UmJ2dzdV58uQJe/DgQYXtU1NTK7xPKzMzk1lZWbFvvvmGRUZGsiNHjjBLS0vm6+sr0+2Rtpr0T0JCAuvSpQvz8vJiMTExLCQkhJmYmLDFixdzdW7fvs2MjY2Zt7c3u3jxIjt9+jRzdnZmDg4OLDc391NuYrWEhYUxPp/PlixZwmJjY9nixYsZn89nZ8+eZYyV/d5v3bolsg2S9OHcuXNZly5d2O7du1lMTAzz9PRkZmZmLDk5+ZNvY23Joo8KCwvZsGHDmLW1NTty5Ai7dOkSmzNnDjM0NGSRkZFy2c6aqkn/fMjDw6PCm4vlsQ9R0pKD7OxsFhQUxCwsLJi5uTmbMGGC2B32Hh4erE+fPhW2ryxpMcbY48ePmZeXF+vatSuztrZm33//fZ3+QK5ITfvn77//ZiNGjGAmJibMzs6OrVmzhhUVFYnV8fDwYObm5szCwoLNmDGDpaWlyXybais0NJQ5OTkxExMT1r9/f3b8+HFumfAD6cqVK1yZJH1YWFjIfvzxR2Ztbc1MTU2Zu7u7wtyEXhFZ9FFmZib7/vvvmb29PTMxMWFDhw5lUVFRn2qTpKq6/fOhypKWPPYhmrmYEEKIwqBrWoQQQhQGJS1CCCEKg5IWIYQQhUFJixBCiMKgpEUIIURhUNIihBCiMChpEakLCgqqcIry8j9BQUHyDrXGjh07BkNDQxw7dkzeoQAANm7cCENDwwqfpyckjPnq1atSe98P52K6evXqJ+uT5ORkBAUFwd7eHiYmJrCxscHEiRNx+fJlmb93TX3K/qnP6IG5ROpGjhwJa2tr7vWNGzdw+PBhjBw5Et27d+fKdXV15RHeZ8vS0hLBwcEwMDCQyvp8fHzQsmVLrFixQirrk9SjR4/g7u4OLS0tuLm5QVtbGxkZGTh16hTGjx+PhQsXwtPT85PGRD4dSlpE6szMzGBmZsa9LikpweHDh9GtWzd88803cozs89a2bVu0bdtWauuLi4vD0KFDpbY+Sa1atQqNGzfGiRMnRJ5Q7uPjg5EjR2L16tUYNGiQwjzQllQPnR4khCiUW7duwdTUVGxKDVVVVYwePRpFRUV49OiRnKIjskZJi8jVxo0b0aVLF0RFRcHW1hZmZmY4evRopddoypdv3LgRZmZmSEhIwLhx49CtWzfY2dlh586dYIzhl19+QZ8+fWBubg4fHx+R9QUFBcHJyQm3bt3CsGHD0LVrV7i6uiI0NFSq25iTk4OlS5fCzs4OJiYm6N+/P3799VexaTEePHiAqVOnwsbGBsbGxrC2tsbs2bORnp4uUu/Zs2eYOnUqLC0tYWVlhXXr1kk0/Xv5a1rC1/Hx8Zg9ezYsLS1hZmaGyZMnV3lt7Pnz5zA0NAQAHD9+XOw6WUFBAZYsWQJra2t069YNXl5eYvOVlZaWIiQkBK6urjAxMYGdnR2WLVsm0Qzb6urquHr1aoXTxA8fPhwPHjwQecp4Xl4e1qxZA1dXV3Tp0gVmZmZwc3PD+fPnxbbp1KlTWLlyJWxsbGBmZobvvvsOWVlZuHv3LkaNGgVTU1O4uLjg3LlzXFvhtarY2FgEBQXB3NwcPXv2xLx585CVlVXlttSmHz5XdHqQyJ1AIMDChQvh4+ODoqIidO/eHWfPnpW4fXFxMby8vNCvXz84OzsjLCwMq1evxpUrV5CWlgYvLy+8efMGu3btwrx587Bv3z6ubXZ2Nnx9fdG7d28MGzYMkZGR+OGHH/D27Vv4+/vXetsKCgrg4eGBly9fwt3dHV9++SWuXLmC5cuXIzk5GYsXLwZQNneTu7s79PT04Ofnh0aNGuHmzZs4efIk/v33Xy7mjIwMjBo1ittmNTU1HDx4EG/fvq1xjJMmTYKBgQFmzpyJ1NRU/Prrr3j16hV+++23CutraWkhODgYgYGBsLCwgJubGwwMDPD06VMAwOrVq2FoaIipU6fi1atX2L17N3x8fBAdHc1Nv75gwQKcOHECQ4cOhbe3N54+fYrQ0FDcvHkToaGhaNiwYaXxDh8+HNu2bcOgQYNgZ2cHe3t79OzZE+3bt4eysrJIXcYY/P398fDhQ3h4eEBXVxfp6ek4dOgQpk6dioiICJFTpqtXr0bLli0xZcoUJCQk4MCBA3jz5g0SExMxbNgwDB48GHv37kVgYCCMjY1FJhX94Ycf8MUXX2DatGl4+fIl9u/fj/v37yMsLKzSuaVq0w+fLZk+jpcQ9t8TpMPCwsSWbdiwgfH5fLZhw4YKy1NTU6ssF75esWIFV+fJkyeMz+czMzMzlpmZyZXPnj2bGRoassLCQsZY2bQKfD6fLVu2jKsjEAjYmDFjWNeuXUWm8ajONpWP19jYmMXHx4uUr1mzhvH5fPbo0SPGGGOLFi1ipqam7M2bNyL1Zs6cyfh8Ple+YsUKZmhoyO7fv8/VycjIYD179qywvyqKWfgkb+HrKVOmiNRbtGgR4/P5LCkpqcpt4/P5bO7cudzrK1euMD6fz4YOHcqKi4u58o0bNzI+n8/+/PNPkXqhoaEi6/vjjz8Yn89ne/bsqfJ9i4uL2dKlS1mnTp0Yn8/nfvr27cu2bNnC/X4ZK5uOpqL3unTpEuPz+SwkJIQx9t/MCfb29uzdu3dcveHDhzM+n88OHjzIlV2+fJnx+Xx2+PBhke3p3bu3yIwKR44cEXlvYT3hPlPbfvhc0elBUif06tWrVu379evH/b9du3YAAHNzc2hpaXHlbdq0AWMMGRkZIm0/PKJSVlbG2LFj8f79e5EJAmsqMjISfD4fLVu2RFZWFvcjjPfChQsAyr6lx8TEoGnTplzbvLw87pt2QUEBAODSpUvo0qULjI2NuXrNmzfHwIEDaxxj//79RV4bGRkBgFg/ScrV1RUqKv+dxBFOSS9cX2RkJHg8Hnr37i3SJ507d0bLli1x8eLFKtevoqKChQsXIjo6GoGBgbC2toaqqiqeP3+O9evXY/To0Vx/mZqa4u+//8awYcO49iUlJSgtLQUAseni7ezsuKNB4L99ycnJiStr06YNAODff/8Vaevu7i4yxfzQoUPRpEkTxMTEVLgdte2HzxWdHiR1QvPmzWvVvkWLFtz/hR+Y5dcpPHUk/MACgKZNm4q0BcCd8klLS6tVTEDZ9af379+L3ALwoZcvXwIAeDwe3rx5g+3bt+Px48d49uwZXrx4wV2rEsaclpYGR0dHsfW0b9++xjGWH2UnPJVVUlJSo/V9+EUBAJcEhNPWP3v2DIwxODg4VNheXV1dovfR0dGBj48PfHx88P79e1y4cAHr16/H/fv3sX//fvj5+QEo2x8OHTqEa9euISUlhfudABC7Flh+n6loX1JSUqqwbYcOHcTatmnTptL9SFr98LmhpEXqBOEHwcdU9kFa/loGUJYIPqZBgwZiZcIEUdE6q6ukpATdu3fHlClTKlzeqlUrAMDFixfx3XffoVWrVujZsyd302xcXBy2b9/O1efxeCgsLBRbT/kP0OqQtO+ltb7S0lKoq6tj06ZNFS6v6jrOjRs3EBkZifHjx0NbW5srV1NTQ//+/WFubo6+ffvi5s2bAIC3b99i1KhRSE1Nha2tLfr27YtOnTpBR0cHI0aMEFv/h0eIH6rpvlRSUlJpf9SmHz5nlLRInST8Qy8qKhIpr+kpq8pkZGQgPz9f5FttcnIyAIhcZK8pHR0d5Ofni4xmA8pGFP7111/ceyxduhR6enoICwvDF198wdU7ffq0SLs2bdpw8X3owydT1HU6OjqIi4uDiYkJGjduLLIsIiJC5BRpeampqdizZw86depU4T1i2tra0NTU5I7u9u7di6dPn2LPnj0iR7vCpCZNz549E3ldXFyMtLS0So+ya9MPnzO6pkXqpJYtWwIA4uPjubK8vDzExsZK9X0YYzhw4AD3WiAQ4Ndff4WmpmalHzbV0bdvX8THx4tdn9i6dSumT5+OJ0+eACgbxdi6dWuRhPXy5UtERkYC+O8I09nZGU+ePMGlS5e4erm5uTh58mStY60uJSUlkVOtkurbty+Asj74UExMDKZNmyaWqMu3VVdXx8aNG8VuBQDKrhO9efOGO4WanZ0NQPTUHWMM+/fvB1D2+5aWw4cPc6dAAeDo0aPIzc2Fs7NzhfVr0w+fMzrSInVSv379sGzZMvzvf/9DWloaVFVVceTIEZEPdWnZsmUL0tLS0LFjR4SHh+PWrVv48ccf0ahRo4+2PX78OG7fvi1WbmRkhNGjR8Pf3x+RkZGYMmUKRo0ahY4dO+LGjRs4efIk7O3tYW9vDwCwt7fHuXPnsGjRInTp0gXPnz/HkSNH8O7dOwD/DRgYN24cTp06halTp8LLywtaWlo4fPhwrU4P1pSWlhauXbuGI0eOVGsgTe/eveHo6IiQkBA8f/4cNjY2SEtLw4EDB9C6dWv4+PhU2rZx48ZYsWIFZs2ahUGDBmHw4MHo1KkTSktLcf36dYSHh6Nfv37cwBR7e3vs27cP/v7++Pbbb1FcXIzw8HDcv38fSkpKYgMxaiM5ORkeHh4YPHgwkpKSEBoaih49elQ6SKY2/fA5o6RF6iQtLS3s3LkTa9aswYYNG9CsWTO4ubmhffv2mDlzplTf65dffsEPP/yA48ePo0OHDti0aZPIaLGqXLt2DdeuXRMrd3R0xOjRo9G0aVMcPnwYGzZswO+//47Dhw+jdevW+O677+Dn58edBhXe4xMTE4OTJ0/iyy+/xJAhQ+Dk5ITRo0fjypUr6Ny5MzQ0NHDw4EGsWrUKhw8fRklJCQYMGICOHTti2bJlUu2XjwkICMCaNWuwdOlSLF26FF999ZVE7Xg8Hn7++Wfs2rULJ06cwIULF6ClpQVnZ2dMnz5dbGBMec7Ozjh27Bh2796NS5cu4dixY1BSUkKHDh2wcOFCjBo1iutXe3t7LFu2DCEhIVixYgWaNGkCY2NjHD58GN9//71UHx4cEBCA27dvY/Xq1dDU1IS3tzemTZtW6TWt2vbD54rH5PEVjZA6ICgoCMePHxd7WgMh1XH16lWMHTsWP/30k8jQeiIbdE2LEEKIwqCkRQghRGFQ0iKEEKIw6JoWIYQQhUFHWoQQQhQGJS1CCCEKg5IWIYQQhUFJixBCiMKgpEUIIURhUNIihBCiMP4PRbqqbjDn6c0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(simulations, bins=20)\n",
    "plt.axvline(0, color='red', linestyle='dashed', linewidth=2)\n",
    "plt.title('Approximate Sampling Distribution')\n",
    "plt.ylabel('# of Simulations')\n",
    "plt.xlabel('Trump Lead in the Sample')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the 100,0000 simulated polls, we find Trump a victor about 60% of the time: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 336,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.60798"
      ]
     },
     "execution_count": 336,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.count_nonzero(np.array(simulations) > 0) / 100000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This number represents the chance that a given sample will correctly predict Trump's victory *even if the sample was collected with absoutely no bias*. In other words, even a non-biased sample will be wrong about 40% of the time. \n",
    "\n",
    "\n",
    "We have just studied the sampling error, and found how our predictions might look if there was no bias in our \n",
    "sampling process. Next we will see what happens when a little biasenters into the mix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "q2",
     "locked": true,
     "schema_version": 2,
     "solution": false
    }
   },
   "source": [
    "## Simulation Study of Selection Bias\n",
    "\n",
    "\n",
    "\"In a perfect world, polls sample from the population of voters, who would state their political preference perfectly clearly and then vote accordingly.\" {cite}`grotenhuis2018`\n",
    "\n",
    "That's the simulation study that we just performed. \n",
    "\n",
    "\n",
    "It's difficult to control for every source of bias. We investigate here the effect of a small, education bias on the polling results. \n",
    "Specifically, we examine the impacts of the 0.5 percent bias in favor of Clinton that we described earlier. \n",
    "\n",
    "This bias essentially means that we see a distorted picture of voter preferences in our polls, where instead of 47.46 percent votes for Clinton, we have 47.96, and we have 48.18 - 0.5 = 47.68 percent for Trump. \n",
    "We adjust the proportions to reflect this bias: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 337,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.48, 0.48, 0.04])"
      ]
     },
     "execution_count": 337,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "proportions_bias = np.array([0.4818 - 0.005, 0.4747 + 0.005, 1 - (0.4818 + 0.4746) ])\n",
    "proportions_bias"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 338,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2939699, 2957579,  268814])"
      ]
     },
     "execution_count": 338,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "votes_bias = np.trunc(N * proportions_bias).astype(int)\n",
    "votes_bias"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, our simulations find Trump winning in about 45% of the samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 339,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulations_bias = [trump_advantage(votes_bias, n) for i in range(100000)] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 340,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Trump Lead in the Sample')"
      ]
     },
     "execution_count": 340,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAEtCAYAAADZSROFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABdG0lEQVR4nO3dd1gUxxsH8O8BInpgQRF/oCCiBwqKKIhIEUEQe4sdEEVBY+81lmgs2MUkllhiQzQK9oCIomjUYO8R6ShGQJCi1Pn9gbfxOMAD7uCE9/M8PsnNze6+O+zdezszu8tjjDEQQgghpEopVHUAhBBCCKGETAghhMgFSsiEEEKIHKCETAghhMgBSsiEEEKIHKCETAghhMgBuUzIf/31FwwMDGBhYYGcnJyqDkcqXF1dYW9vX+nbTU5ORlZWltTWV1BQAF9fXwwbNgwdO3ZEhw4d0KdPH2zevBnp6elS2460FW3/BQsWwMDAoFJjuHXrFgwMDET+tW3bFp07d8bo0aNx6tQpsWVOnjwJAwMD3Lp1q8zbi4uLk6iegYEBFixYUOJraSgaS1V9Hvbv3w9ra2u0b98eGzZsqPTtf42Pj4/YMWJkZAQbGxvMnz8fiYmJxdaPj4+voohLVpZj99GjR5g2bRqsrKxgbGwMa2trzJw5E48ePaqESMunIp/NkihJbU1SdPbsWdStWxepqakICQmBs7NzVYdUYRMnTsTHjx8rdZuhoaGYM2cO/P39UbduXamsc968eTh//jx69eqFfv36QUFBAY8fP8Zvv/2GP//8E0ePHkXDhg2lsi1ZGj58OCwtLatk246OjnB0dAQA5OXlITk5GcHBwZg3bx7u3r2LFStWcHXNzc3h7e0NfX39Mm1j6dKliIqKwsGDB79a19vbGzo6OmXbiTI4ceIEVqxYgYcPH3JlVfF5ePHiBdasWYMOHTpg+vTpMDQ0rNTtl8XEiRPRsmVLAEBOTg7i4+Nx9OhR3LlzBwEBAVBVVQVQeCzp6OhAXV29KsOtkKtXr2LixIlo1aoV3NzcoK6ujjdv3uDEiRP4888/4ePjgx49elR1mJVC7hJyTk4OgoKCMGDAAJw9exb+/v7VIiFbWVlV+jYfPnyIDx8+SG19d+/exZkzZ7BgwQKMHTtW5D1bW1vMmDEDv/32G+bOnSu1bcqKqakpTE1Nq2TbBgYGGDBggEjZ+PHjMX/+fBw9ehQWFhbo3bs3AKB58+Zo3rx5mbcRFhYGbW1tieoWjUXa/v77b2RnZ4uUVcXn4Z9//gEAeHl5VcnZeVl07doVFhYWImUdO3aEp6cnAgIC4OLiAgAwNDSU6x8Wkvjpp59gaGgIPz8/1KpViyt3c3PDgAED8OOPP8LOzg5KSnKXrqRO7rqsQ0ND8eHDB1hYWMDa2hrXrl3Du3fvqjosAuDevXsAiv8y7dWrFzQ1NXH//v1Kjqp6UFBQwLJly1C/fn3s3r27qsOplnJzcwEAfD6/iiMpny5dugAAIiIiqjgS6UlJSUF0dDQsLCxEkjEANGjQAAMHDkRSUpJcdsnLgtwl5DNnzoDH48Hc3ByOjo7Iz88XG1uzt7fH4sWLcfz4cTg4OKBDhw4YMWIEbt68We56S5YswaJFi9CuXTvY2toiJSUFABAeHg53d3fujMrNzQ1///03t+zVq1dhYGCA6dOni6zzhx9+gIGBAa5evQpAfMzM1dUVXl5eCA4ORv/+/dGuXTv06dMHoaGhyMjIwNKlS9G5c2dYWlpi6dKl+PTpE7csYwy+vr747rvvYGpqinbt2sHZ2Rm7du2C8E6oCxYswPbt2wEADg4OcHV15ZaPiIjA5MmTYWZmBhMTE4wYMQLXrl376t9G+EV27NgxFBQUiL0fHByMw4cPi5T99ddfGD9+PCwsLLixsKVLl4qcuS9YsAB9+/bFnTt3MHz4cLRv3x4ODg7w9/dHbm4uNm7cCCsrK3Tu3BkzZszA+/fvRdrR3d0dISEh6N27N9q3b4+BAwciMDCw1H0pOoa8YMECODs74+HDh3BxcYGJiQm6du2KVatWibQ9AERGRmLSpEkwMzODhYUFVq1ahWPHjlV4LE9VVRXdu3fH06dPkZSUBKD4carAwEAMGTIEpqam6NSpE8aOHYs7d+5w7xsYGCAhIQG3b9+GgYEBTp48ifj4eBgYGGD//v0YOXIkjI2N4e7uztUvbsx4x44dsLGxgYmJCdzc3ES6nEtb7styV1dX+Pv7F1te9Cz1xYsX+P7772FmZob27dtj2LBhCA4OFqnj6uoKDw8PXL16FYMHD0a7du1gZ2cHHx+fYo/JL5dbuHAhgMIzry//9mXZ7ubNm2FqagpLS0u8ePGixO39+eefcHFxQadOnWBsbAx7e3t4e3tXaE7MmzdvAECkx6S4MeQnT55g6tSp6Nq1K4yMjGBpaYnZs2eLjT/7+vqiX79+MDExgYWFBSZPnoyXL1+K1MnOzsbmzZthb28PY2NjODg4YOvWrWL7kZycjIULF6JLly7o1KkTli5dKtG+1qlTB4qKirh06VKxJ15Tp07FkydP0KJFC67s3bt3+PHHH+Hg4ABjY2N06tQJbm5uIp8B4XyNGzduYPHixTA3N0enTp2wcOFCZGVlITQ0FAMGDICJiQkGDBiAv/76i1tW+Jl7+PAhJk2ahA4dOsDa2hpr1qwR+y4oStL2Kolc9QFkZGTgypUr6NChAxo3boxu3bpBWVkZ/v7+GD9+vEjdGzdu4PTp03B1dYWGhgZ8fX0xfvx47N27F507dy5zvXPnzkFPTw+LFy9GUlIS1NXVcenSJUyZMgU6OjqYNGkSAOD48eNwd3fHtm3b4ODgAFtbWwwaNAj+/v64du0abGxscP36dRw7dgwjRoyAra1tifv75MkT3Lt3D25ublBTU8POnTsxY8YMtGnTBnXq1MHMmTMRHh4OPz8/NGnSBFOmTAEAbNmyBTt27MCgQYMwbNgwZGZmIiAgABs3boSGhgYGDRqE4cOHIyMjAxcvXsTChQvRunVrAIVfPqNGjULjxo3h5eWFWrVq4ezZs/D09MTGjRu5rtLiODk5YdOmTTh48CBCQkLQs2dPWFpawszMDHXr1oWysrJI/bCwMEyYMAEdO3bEtGnTwOPxcP36dfj5+SE3Nxdr1qzh6r579w4TJ07E0KFD0b9/fxw4cACLFi3CmTNnkJ6eju+//x6vXr3C4cOHUadOHZFlX716hWnTpmHIkCEYMWIEAgICMG3aNGzYsAH9+vUrcX+KSklJgYeHB3r16oX+/fvj6tWrOHjwIJSVlTFv3jwAwOvXrzFq1CgAwLhx46CkpITDhw/jzJkzEm+nNF/+nRo3biz2/u3btzFz5kzY2tpi6NCh+PjxIw4dOoSxY8fi3LlzaN68Oby9vbFmzRo0bNgQEydORMeOHbnlt27dim7duqFfv36oXbt2iXEEBgZCWVkZbm5uqFWrFg4cOAA3NzccP36ci1ESEydOREFBAcLDw0sdq3748CHc3NygqqqKsWPHgs/n49SpU5g8eTKWLl2K0aNHc3X/+ecfzJgxA8OHD8fw4cNx9uxZbN++Herq6iL1isahp6cHPz8/kfHZsmz37t27iImJwdy5cxEfH49WrVoVu63jx49jyZIlsLe3x5w5c5Cbm4uLFy9iz549qFu3Lvc5Lk16ejp3UpCfn4+4uDh4e3tDS0sLQ4YMKXE54edbV1cXnp6eqFOnDu7evYtTp07h33//5eYUnD59GsuXL8fAgQPh6uqKlJQU/P7773B1dcXFixehpqaG/Px8eHl54e7duxg2bBj09fXx+PFj7NixA8+ePcOvv/4KHo+H7OxsuLi4ID4+Hm5ubtDQ0IC/vz/Onz//1f2sU6cOevfujTNnzqBHjx6wt7eHtbU1unTpAm1tbbFu6k+fPmH06NFIT0/H6NGjoampiejoaPj6+sLT0xOhoaHc+DpQ+EO7VatWmD17Nm7fvo2TJ08iMTERT58+haurK9TU1LBr1y5Mnz4dwcHBqFevHrfs9OnT0aRJE8yePRvPnj3D/v37ERERgT179hS7L5K2V6mYHPnjjz+YQCBge/bs4co8PT2ZQCBgDx484Mq6d+/OBAIBu3jxIleWnJzMzMzM2LBhw8pVz9DQkMXExHBlubm5zNbWlnXr1o2lp6dz5WlpaczGxobZ2NiwnJwcxhhjqampzMrKijk6OrKUlBRmZ2fHevTowTIzM7nlXFxcWPfu3UVeCwQCFhISwpUdOnSICQQCkdgKCgqYra0tGz58OGOMsZycHNaxY0c2c+ZMkbZLT09nxsbGzMvLiyvbtm0bEwgELC4uTmS7RWPLzc1lo0aNYl27dmXZ2dmsNPfv32cODg5MIBBw/4yMjJiXl5fI34gxxjw8PFj37t3F1jls2DBmamrKvZ4/fz4TCATs4MGDXNmVK1eYQCAQW37EiBHM2tparB337dvHlX38+JE5Ojoya2trlp+fz9X7sv2F2yz6+sCBAyKx9urVS2R7CxcuZG3btmURERFcWWJiIuvQoYNYWxd18+ZNJhAI2LZt20qsc+zYMSYQCNjZs2cZY4ydOHGCCQQCdvPmTcYYY8uWLWOmpqasoKCAW+b58+fMycmJXbhwgSvr3r07c3Fx4V7HxcUxgUDAHB0dueNWSCAQsPnz54u8btOmDXv+/DlXFh0dzYyMjNiUKVNKXK6k8qJtzZj432Po0KGsQ4cO7M2bN1zZp0+f2KBBg1j79u1ZcnIyt5xAIGCXLl0SqWdubs59RkpStC3Ls90vly2Js7MzGz58uMjfSPh90rdv31KXFX5mi/tnaGgo8l32ZX3hcbd06VJmYmLC3r9/L1Jv5syZTCAQcOXjx49nffr0Ealz5coV1rt3bxYeHi7SXlevXhWpd/ToUZHv1YMHD4p9z2ZmZrLevXtL1GYZGRlsxowZYvvbp08fdujQIe4zzBhj586dKzYmX19fJhAIWGBgIGPsv8/akCFDuOXz8/OZlZUVEwgELDQ0lFtW+JkLCwsT2e8hQ4aIfPds2rRJZNtFjydJ26s0ctVlLTzLEM5A/fL/hd1eQi1bthSZeaeuro4BAwbgwYMHSE5OLnM9HR0dkV/vT58+RWJiIkaPHi3yi6tevXpwcXHB27dv8fjxYwBA/fr1sWLFCsTExGDo0KFITEzEunXrvjqzuXbt2rCxseFe6+npASjsYhbi8XjQ1tbmunNq1aqFGzdu4McffxRZ1/v376GqqlrqJU7v37/H7du30a1bN3z69AkpKSlISUnBhw8f4OjoiKSkpK9eZmBiYoI///wTO3fuxPDhw9GsWTPk5ubi8uXLGD58uMiZ4s6dO3HixAmRM+fS4vzy7y7sorKxsRFZvlmzZmJdW2pqatxZKwCoqKhg5MiR+Pfff7m/kaR69eol8trQ0JA7ThhjuHTpEmxsbERmPWtqaqJ///5l2k5JhOOcJf2Sbtq0KTIzM7Fq1Sq8evUKQGFXcGBgoESTH7t06SI2VlccGxsbkW5dXV1d2NjYICwsDPn5+ZLsisSSkpLw4MEDDBgwAE2bNuXKa9euDQ8PD3z69Ak3btzgyuvUqQM7OzuRenp6elw3v6y2q6KiAnNz86+u9/Tp09i1a5fI3zA5ORn16tWT+BLE+fPnY9++fdi3bx92796NVatWwcTEBFOmTBH7LvzS8uXLERISggYNGnBlGRkZXG+IcPtNmzZFZGQktm/fznV3d+vWDefOnUOnTp0AAEFBQVBXV4eRkRH3XZGSkoJu3bpBUVERV65cAVA4bNe4cWOR79m6deti6NChEu0rn8/H5s2bcf78eUydOhWmpqZQUlLCy5cv8eOPP+L777/njrnevXvjr7/+grW1Nbf8l93BRdvXwcEBCgqFaU5BQQHNmzeHioqKSM9ls2bNAEDse2XcuHEi3z3CiawhISHF7oek7VUauemy/vfff3H79m20aNECPB6PO0gMDQ3B4/Fw7tw5LFy4kGug4rqLdHV1wRhDQkICGjVqVKZ6wv8KCbcvTJJfEnZ3vX79mpup6+DgACcnJwQFBWHkyJEi3YQladCggUiXjKKiYrGxKCoqcmPDQGFSvnLlCi5duoSoqCjExMQgLS0NAETqFSW8FvTgwYMlXg4jHKcqjZKSEuzs7LgvxcjISBw5cgQHDx7EqlWr4OjoCBUVFSgqKiIuLg5bt25FREQEYmNj8fbt2xLX++V+S9oWQOGPqaLd5bq6ugCAhIQEtG/f/qv7JFT08hFlZWXuyyA1NRWpqaki41lCwmOiolJTUwGgxEvHXFxcEBYWhkOHDuHQoUNo1qwZunfvju+++06i2baSXh5T3P7o6OggJCQEKSkp0NDQkGg9kkhISABQ/GdN+MPn9evXXFmDBg24L1khZWXlUseQZbXd4tSqVQt///03zp49i8jISMTGxnI/6iSd+W5kZCQ2y7p///7o168f1q5di169ekFFRUVsOR6Ph/fv32Pnzp148eIFYmNj8fr1a+4zI2yjyZMn4/79+/Dx8YGPjw9atWoFe3t7DB06lDsxiY2NRUpKSomXBwq/KxISEoq9EqC4di2Nvr4+pkyZgilTpiA9PR1//vkntm7disuXLyMwMJAbTuPxeNi1axfu3buH2NhYxMbGcj9kix4DRYd9lJSUxD4Dwr9p0WWLXmrYoEEDNGjQgDtuipK0vUojNwn5/PnzyM/PR3R0tMgZolBaWhqCg4O5P0pxv/KFX5zCL/Oy1Pvy/4HSE5vwvS/XnZWVhadPnwIArl27hqysrK+eIZc0jb+0cQbGGObOnYuzZ8+iU6dOMDU1xfDhw2Fubo4xY8aUuj3hfo8ePbrE6/pKGhcDgO3bt0NTU1Psl2/Lli2xZMkS5Obm4ujRo4iIiICxsTGOHj2KZcuWQU9PD2ZmZnBycoKJiQkOHjxY7Jhrce3x1TEXFP83Fn64iv5dv6a0L9y8vDwAEEv+AEodjy2LZ8+egcfjlXjTElVVVRw6dAj3799HcHAwN859+PBheHt7f3XMvKzt8SVJ2rQ8Z8+lfdaE2/zybyxJUpTFdiVtu40bN2LXrl1o27YtOnTogAEDBsDU1BQrV66U6Eu5JLVr10b37t2xf/9+REZGom3btmJ1rly5gu+//x5NmjRBly5dYGtrC2NjY4SFhWHnzp1cvaZNm+LUqVO4desWLl26hGvXrmHXrl3Yt28fN78mPz8fLVq0wLJly4qNRzjeKhxHLqq09v0y3uvXr2Pu3Lkinys1NTUMHToUAoEAw4YNw507d9C7d28kJCRg+PDhyMrKgrW1NXr37o02bdqAMYbJkyeLrb+4v5kk3ylAybmjpONP0vYqjdwkZOHs6rVr14p0EQPA8+fP4ePjA39/fy4hx8bGiq0jJiYGioqKXBdEWeoVJfwlGxkZKfZeVFQUAIh0c23atAkJCQmYN28e1q9fj02bNmHJkiWl7XK5hIeH4+zZs/j+++9FZnbn5eUhNTW11GtWhfukqKiIrl27irwXERGB+Ph41KlTp8TlAwICAADfffddsQe1QCAAUNilmJ2djbVr18LCwgJ79+4VSbZbt279+o6WQXx8PBhjIjFFR0cD+O9MWRoaNWqEunXrcuv+UkxMTIXXn5GRgbCwMJiampZ4JhsVFYX09HR06NABHTp0wJw5cxAREYHRo0dj3759ZZrEVprizgJiYmKgpqbGnb0rKCiIzR4ta7cxUPbPmrTIYrsJCQnYtWsXBgwYAG9vb5H3ytM2RQl/KJSUFFauXAldXV2cOHFC5ISg6A9g4QxxS0tL7ozuzp07GDNmDA4ePIjOnTujWbNmePz4Mbp06SKyPeEkNWHbNGvWDOHh4cjLyxP5nEtyp7gnT57gwIEDcHR0FJlkKyScQCjsDdi+fTuSk5Nx4cIFkZ4qaU2q/FJcXJxIT1FKSgrS09OL7SEDIHF7lUYuxpCjo6Px+PFjdO7cGQMHDkSPHj1E/nl5eUFDQwPXr1/nujwfPXokcs1rUlISTp8+jS5duqB+/fpcuaT1ijIyMuJmZWdkZHDlGRkZOHLkCDQ0NGBsbAyg8EA+fPgwhg0bBg8PDwwZMgSHDh1CeHi4lFroP8IuzaJnsseOHcPHjx+5szjgvw+t8JdqkyZNYGxsDH9/f5Gu49zcXCxatAjTpk0TWb6ofv36IS4uDjt27BB7Lzs7GwEBAWjRogVatmyJT58+4ePHj2jRooXIh/TZs2e4ffs2AJS6rbJISkrChQsXuNcfP36Er68vWrRoIdXbYyooKMDe3h5Xr14V+bJJS0vD2bNnK7RuxhhWr16NrKwssSsKvrRq1Sp8//33yMzM5MpatmyJevXqiXwJKCgolLkL90vXrl0TOUb++ecfhIWFwd7envvh07hxYzx//lzkTKi4mbUldQkKCT9Lp0+fFrk0JycnB/v27YOysrJMbiQii+0Kh46Kfj5DQ0MRHR1doWP+06dPuHTpEtTV1UvsyUpNTYWWlpZIMn7z5g2CgoIA/NeDMX36dMybN0+kR6Nt27aoVasW9/eyt7dHamoqfH19RbZx9OhRzJw5k7tUyMnJCenp6Th+/DhXJzc3F8eOHfvqPvXp0wcKCgpYt25dsTcx8vPzA/DfvJrU1FTUqVMHWlpaXJ2cnBwcPXpUZP+k4dChQyLHtnB29ZdzXb4kaXuVRi7OkIW/br777rti369VqxaGDBmCHTt2cNckKysrY8KECRgzZgxUVFRw5MgRFBQUcJenCElar7ht/vDDD5gxYwaGDBnCxfbHH3/g33//xbZt26CgoIDs7GwsXrwY6urqmDNnDgBgzpw5CA4OxuLFi3Hq1Klix3rKy9TUFKqqqlizZg1ev36NevXq4datWzh//jxq164t8kUtPMv67bffYGtrCwcHByxZsgRjxozBkCFDMHLkSDRo0ADnzp3DgwcPMHv27FJve+nl5YVbt25hy5YtCA0NhYODA3ebuzNnziAxMRF79+4Fj8dD/fr1YWJigpMnT0JVVRV6enp4+fIljh8/zn3gMzMzS/1RJKlatWph4cKFePLkCZo0aYITJ07g7du3xf5wqKjp06cjNDQUw4cPh6urK5SVlXH06FHuy0SS7rAXL15wx3F+fj6SkpIQHByMBw8ewM3NrdghG6GxY8diwoQJGD16NAYOHIjatWsjODgYsbGxWLduHVdPXV0dz58/x5EjR9C5c+cyH4PKysoYNWoUXF1d8fHjR+zfvx/16tXDjBkzuDp9+/bF3r17MWXKFNjZ2eHJkye4cOGC2Nm98PW2bdtgYWFR7Bib8Lj87rvvMHLkSPD5fJw+fRpPnjzBkiVLJOruKw9pb7dVq1bQ0tLCjh07kJ2djaZNm+Lhw4fw9/cX+3yW5saNGyI/ElJSUnDixAkkJCTgxx9/LHG4y9bWFufPn8fSpUvRrl07xMfHcz/WAXDb9/DwwJIlS+Du7g5nZ2cwxnDq1ClkZ2dzEySHDh0Kf39/rFy5Ek+ePEH79u3xzz//wM/PD0ZGRhg8eDCAwju9HTt2DCtXrsSrV6/QokULnD59WqIbOrVo0QILFy7E6tWrucsNhT/or1+/jsuXL8PV1ZWbk2Nra4uQkBB4eXnB2dkZ6enpCAgI4HpCJW1fSdy6dQsTJkxA9+7d8eDBA5w6dQoDBw7kJr0VJWl7lUYuEvLZs2ehpqYGJyenEusMGzYMu3bt4mYYCh9q8MsvvyA9PR1mZmaYPXu22MQWSesVp2fPnti7dy9++eUX/Pzzz1BSUoKJiQl++uknmJmZASi8MD8qKgrr16/nPrwNGzbE3LlzsXjxYmzduhXz588vb9OIady4MXbt2oUNGzbgl19+gbKyMvT09LBp0yY8fPgQBw4cQFJSEho3bow+ffogKCgIJ0+exO3bt+Hg4ABTU1P4+vrCx8cH+/btQ15eHvT09LB27VoMGjSo1G2rqKjgwIED8PX1xYULF/Dbb78hMzMT6urq6Nq1K7y8vEQmcmzduhVr1qzBiRMnkJOTA21tbXh6ekJfXx9Tp07FzZs30bNnzwq3SZMmTbBo0SKsW7cO7969g5GREfbt2yfRjNiy0tHRwaFDh7Bu3Trs3LkTtWvXxsCBA6GoqIg9e/YUO75c1MWLF3Hx4kUAhT8mmjRpAj09PWzevLnU68ABwNraGr/++it27tyJX375BdnZ2WjdujU2bdqEPn36cPWmTp2KZcuWYfXq1Zg8eXKZu7KHDx8OHo/HJRYLCwssWLBA5Mxk+vTpyMvLw7lz5xAWFgYTExP8/vvv3A9ToZEjR+LmzZv47bff8OjRo2ITsvC43LZtG/bu3YuCggIYGhri559/lul9jKW9XWVlZezatQtr167FgQMHwBiDjo4OFi1ahLy8PPz00094/Pgx17tWki9/TCooKEBNTQ2GhobYunVrqbPply9fjrp16yIkJASnTp1C06ZNMXDgQDg6OnJ/h7Zt22Lo0KHc9eWbNm1CQUEBjI2NsXv3bm4ymbKyMvbv34+ff/4ZgYGBOH36NJo0aYKRI0di8uTJ3PCWoqIifvvtN2zevBkXLlxAVlYWbG1t4e7ujpkzZ361zdzc3NC2bVscPnwY58+fR0pKClRUVGBoaCh2XI8YMQIfPnzA8ePHsWrVKjRu3BgdOnTA9u3buZs+CW94U1GrV6+Gv78/1q1bBw0NDcyePbvU3itJ26s0PCbJyLucsbe3h7a29ldvnC9pPfLtcnV1RUJCQomXIkhbcnIy1NXVxc6EV65cCV9fXzx48ECiy4oIIfLp5MmTWLhwIQ4cOCA2013W5GIMmZBvxfTp09GnTx+R8dCPHz/i8uXLMDQ0pGRMCCk3ueiyJuRbMWDAACxZsgSenp5wcHBAdnY2Nynoy8cmEkJIWVFCJqQMhg4ditq1a+PAgQNYv349FBQUYGxsjP379xd72QYhhEjqmxxDJoQQQqobGkMmhBBC5AB1WcsQYwyy7H8QTvSlPo7iUfuUTuFD4U0sCupV/Frw6oaOndLJsn14PMlvb1ndUEKWIcaA5OSMr1csp/r1C69rS0v7KLNtfMuofUqn0aTwJjDJ/4rfIammo2OndLJsn0aNVFFD8zF1WRNCCCHygBIyIYQQIgcoIRNCCCFygBIyITVUTnYucrJzqzoMQshnlJAJIYQQOUAJmRBCCJEDlJAJqaGUunSGUhe63Sch8oKuQybkG8Xn14aSUvl/UyvcuyfFaAghFUUJmZBvlJKSAnLyGaIS0sq8rJ52fSjLICZCSPlRQibkGxaVkIZFv14v83KrJ1mhnQziIYSUH40hE0IIIXKAEjIhhBAiByghE0IIIXKAxpAJqaHY+PEoKKDnCxIiLyghE1JDsZ27kJ+XD9AjBgmRC9RlTQghhMgBOkMmpKa6cwe8/HxAv21VR0IIASVkQmoshc7mhV1k/36o6lAIIaAua0IIIUQuUEImhBBC5AAlZEIIIUQOUEImhBBC5AAlZEIIIUQOUEImhBBC5IDcJORnz57ByMgIiYmJIuWOjo4wMDAQ+5eSksLVefToEVxdXWFqagpra2ts2rQJubm5IuuJjo7GxIkTYWZmBgsLCyxbtgwZGRkidZKSkjB79mxYWFigU6dOmDVrFt69eye7nSakChXc/hu5f92s6jAIIZ/JxXXIkZGR8PLyQl5enkh5ZmYm4uLiMHv2bHTu3FnkvXr16gEAYmJi4O7uDlNTU2zZsgWvXr3C5s2bkZGRgaVLlwIA0tLSMGbMGGhoaGDdunVITk7G+vXrkZiYiJ07dwIA8vLy4OHhgaysLCxfvhx5eXnYuHEjxo8fjxMnTkBJSS6aihDp6dQJjG6dSYjcqNIsk5eXBz8/P2zcuBG1atUSe//FixdgjMHBwQH6+vrFrmPXrl1QU1PDL7/8AmVlZXTr1g0qKipYtWoVvLy8oKmpicOHD+PDhw8ICAhAw4YNAQCamprw9PTEgwcPYGJignPnzuH58+c4f/48t602bdqgb9++CAoKQu/evWXXEIQQQmq8Ku2yvnPnDjZs2IBx48Zhzpw5Yu8/e/YMtWvXRosWLUpcx/Xr19G9e3coKytzZc7OzsjPz0dYWBhXx9zcnEvGAGBtbQ0+n4/Q0FCuTqtWrUQSv/C1sA4h1cX/GvPB8/JErcmTUL9+nTL/4/NrV/UuEFLtVOkZsr6+PoKDg9GoUSOcPHlS7P0XL16gQYMGmDVrFq5fv478/HzY2dlh0aJF0NDQwMePH/HmzRvo6emJLKeurg5VVVVERUUBKOwS79+/v0gdRUVFNGvWTKRO0fUAgI6ODlenrHg8oH79OuVaVhJKSooAZLuNb1l1bx/h/pWHSm0l8H77DQDwYv6aMi2rp10fdVVqQUlJbqagSF11P3YqSpbtw+NJfZXfjCpNyI0bNy71/efPnyMpKQmtW7eGq6srIiMjsW3bNri5ucHf3x/p6ekAAFVVVbFl+Xw+N2krPT1dojqtWrUqtk5MTEyZ942Qb8WiX6+Xqf7qSVYwatlIRtEQUnPJ9UylJUuWgDEGExMTAICZmRn09fUxatQonD59Gt26dQMA8Ir5ScUYg4LCf7/gpVWnLBgD0mQ4YUb461SW2/iWVff2qcqzt7y8/GrbrkD1P3YqSpbt06iRao09S5brPqf27dtzyVioU6dOUFNTw/Pnz7mz3qKXLwFAVlYW1NTUABSeQRdXJzMzk1uHJHUIIYQQWZHbhJyVlYUTJ07g+fPnIuWMMeTm5qJhw4bg8/nQ1NQU61JOTk5GRkYGNyasp6cnVic/Px/x8fGl1gGA2NjYYseWCSGEEGmS24Rcu3ZtrFu3Dtu3bxcpv3TpEj59+sRdl2xlZYXLly8jJyeHqxMYGAhFRUWROrdu3UJqaipXJywsDFlZWejatSuAwlnXL1++RGRkJFcnIiICkZGRXB1CCCFEVuQ2ISsqKmLSpEm4ePEiVq1ahRs3bmD//v2YP38+HBwcYGFhAQAYP3483r17B09PT1y+fBn79u3DmjVrMGzYMGhpaQEARo0aBWVlZbi7u+PixYs4fvw45s6dC1tbW3Ts2BEA0Lt3b+jq6mL8+PE4d+4czp49iwkTJqB169bo1atXlbUDIbLy0ag9Ipq0rOowCCGfyfWkrrFjx0JVVRUHDhzA8ePHUb9+fYwYMQJTp07l6ujr62Pv3r3w9vbGtGnT0LBhQ4wdO1akjrq6Og4cOIDVq1djzpw54PP5cHZ2xrx587g6ysrK2LdvH3766ScsWbIEysrKsLKywoIFC+guXaRaigi4VOYZ1oQQ2eExxlhVB1FdFRQwJCeLTxSTFpoJWjp5bx8+v3aFruVVUlLEk8jkciVV31W9EZWQVq5lV0+ygoFOA7ltV2mQ92Onqsl6lrWCQs2cZk2nfoRUESUlBeTkM0QlpJVr+bZ66lKOiBBSlSghE1KFynuWChSe5VZEu9YaOAOg36yACq2HECIdcjupixBCCKlJKCETQgghcoASMiGEECIHKCETQgghcoASMiGEECIHypSQMzIycO/ePe51eHg4pk2bhpkzZyI8PFzqwRFCCCE1hcSXPUVERMDNzQ2NGjXCmTNnEBcXh7Fjx4Ixhlq1auHixYvYvXs3LC0tZRkvIURK4lduREBoRFWHQQj5TOIz5C1btgAA5s6dCwA4fvw48vLycPDgQdy4cQNt2rTBr7/+KpMgCSHS936EGwLb96zqMAghn0mckP/++2+4u7vD1tYWABASEgJdXV2YmpqiTp06GDhwIB4/fiyzQAkhhJDqTOKEnJ2djYYNGwIAEhISEBERARsbG5E6ioqK0o2OECIzDY8eQM+HgVUdBiHkM4kTso6ODu7evQsA8Pf3B4/Hg4ODAwCAMYY///wTurq6somSECJ1zX6YjSnBNMxEiLyQeFLXyJEjsWLFCjx+/BiRkZFo3bo1unTpgn/++Qfz58/H8+fPsXbtWlnGSgghhFRbZUrIfD4fZ8+ehampKSZPnsy99+nTJ6xcuRIDBgyQSZCEEEJIdVempz31798f/fv3FykTCAS4cOGCVIMihMiv/zXmQ0lJkXsmblnl5RUgMzNbylER8u0r8+MXs7OzkZqaivz8/GLf19LSqnBQhBD5pVJbCVnZeeV6jrOedn0oK9ENAgkpjsQJOTU1FStWrMDFixdLTMYA8OzZM6kERgiRX+V9jvPqSVYw0Gkg/YAIqQYkTshr167FhQsXYGNjgzZt2kBZWVmWcRFCCCE1isQJOSQkBEOHDsXKlStlGQ8hpJI8evmuXGe5hBDZkHgwJy8vD+3atZNlLIQQQkiNJXFCNjc3x61bt2QZCyGEEFJjSZyQFy1ahPDwcHh7e+Phw4eIj4/H69evxf4RQr4NrQY6YPOhWVUdBiHkM4nHkPv164eCggLs3bsX+/btK7EezbIm5NtQ58lDtKrqIAghHIkT8oQJE8Dj8WQZCyGEEFJjSZyQp06dKss4CCGEkBqtzHfqunbtGoKDg/H69WvUqlULWlpasLOzg7W1tSziI4QQQmoEiRNyQUEB5syZgwsXLoAxhnr16qGgoAAZGRk4fPgwnJycsGXLFurWJoQQQspB4lnWv/32G86fP4+RI0ciLCwMt2/fRnh4OMLCwuDi4oLAwED8/vvvsoyVEEIIqbYkTsgnT55Ejx49sHTpUjRu3Jgrb9y4MRYvXgxHR0f88ccfMgmSECJ9KcNd8Wc7x6oOgxDymcQJOSEhAVZWViW+b2lpibi4OKkERQiRvYRVm/Cz4+SvVySEVAqJE3LDhg0RHR1d4vvR0dFQU1OTRkyEEEJIjSNxQra3t4evry9CQkLE3rt06RKOHj0Ke3t7qQZHCJEdlccPoP82oqrDIIR8JvEs6xkzZuCvv/7C5MmToa+vDz09PQBAZGQkIiMjoa2tjRkzZsgqTkKIlLUe1ANbAPSbFVDFkRBCgDKcITdo0ADHjx+Hh4cHGGO4evUqQkNDUVBQgLFjx+LEiRNQV1eXZayEEEJItVWmG4PUq1cPc+bMwZw5c2QVDyGEEFIjlZiQX79+DXV1daioqHCvJaGlpSWdyAghhJAapMSE7ODgAG9vb/Tr1w9A4aQuSe7CRU97IoQQQsquxIQ8efJkGBgYiLym22ISQgghslFiQp4yZYrIa0me9pSTk1PxiAghhJAaSOJZ1g4ODrh06VKJ7589exY2NjZSCYoQInsv/YMxY/SGqg6DEPJZiWfIKSkpePXqFfc6ISEBjx49Qr169cTqFhQU4OLFi3SGTMg35JOxCV5dy6jqMAghn5WYkGvXro3Zs2fj3bt3AAAej4edO3di586dxdZnjKF3796yiZIQQgip5kpMyHw+H7/++iv++ecfMMawaNEiDBs2DKampmJ1FRQUoK6uDktLS5kGS4i84fNrQ0lJ4pEfEUpKilKOpmy0l8zC5KeJ9IAJQuREqTcGMTIygpGREYDC65CdnJwgEAgqJTBCvgVKSgrIyWeISkgr87Jt9ar2znbqfgfhDFBCJkROSHynrqKzrovz9OlTtG3btkIBEfKtiUpIw6Jfr5d5Od9VNMRDCPmPxAk5NzcXu3btQlBQELKyslBQUMC9l5+fj8zMTGRkZNCNQQghhJBykHjwa8uWLfDx8UFaWhrq1KmDhIQE/O9//4OSkhISExORm5uLxYsXyzJWQgghpNqSOCH/+eef6Ny5M0JCQrB7924AwNKlSxEYGIidO3ciLy8PtWrVklmghBBCSHUmcUJ++/YtnJycoKCgAE1NTTRq1Aj37t0DAHTr1g2DBg3CsWPHyh3Is2fPYGRkhMTERJHysLAwDBkyBCYmJrC3t8fevXvFln306BFcXV1hamoKa2trbNq0Cbm5uSJ1oqOjMXHiRJiZmcHCwgLLli1DRoboNZhJSUmYPXs2LCws0KlTJ8yaNYu77IsQQgiRJYkTsoqKisgZsI6ODv755x/udfv27REXF1euICIjI+Hl5YW8vDyR8rt372LixIlo2bIlfHx80K9fP3h7e2PPnj1cnZiYGLi7u6N27drYsmULxo0bh3379mHNmjVcnbS0NIwZMwZJSUlYt24dZs+ejfPnz2P27Nlcnby8PHh4eODhw4dYvnw5li9fjrt372L8+PFicRFSHXw0ao+IJi2rOgxCyGcST+pq06YNrl69iuHDhwMAWrZsyZ0hA4Vn0GV9+EReXh78/PywcePGYru7t23bhrZt22L9+vUAAFtbW+Tl5WHHjh1wdXWFsrIydu3aBTU1Nfzyyy9QVlZGt27doKKiglWrVsHLywuampo4fPgwPnz4gICAADRs2BAAoKmpCU9PTzx48AAmJiY4d+4cnj9/jvPnz0NfX5/b5759+yIoKIhuekKqnYiAS+WaHU4IkQ2Jz5BHjx6NS5cuYdSoUcjIyECfPn3w9OlTLFy4ELt378b+/fvRrl27Mm38zp072LBhA8aNG4c5c+aIvJednY3w8HA4OTmJlPfs2RMfPnzA3bt3AQDXr19H9+7doayszNVxdnZGfn4+wsLCuDrm5uZcMgYAa2tr8Pl8hIaGcnVatWrFJWMA3GthHUIIIURWJD5DdnZ2xsqVK7Fv3z7UqVMHXbt2xYQJE7gJXlpaWli4cGGZNq6vr4/g4GA0atQIJ0+eFHkvLi4Oubm50NPTEynX1dUFAERFRcHExARv3rwRq6Ourg5VVVVERUUBKOwS79+/v0gdRUVFNGvWTKRO0fUAhV3zwjplxeMB9evXKdeykhDe6UmW2/iWVUb7VPXdtr5FSkqKcn/M0merdLJsn5r8lF+JEzIADB06FEOHDuVez549GyNHjkRaWhr09fVFzlIl0bhx4xLfS09PBwCoqqqKlPP5fABARkZGiXWE9YSTttLT0yWq06pVq2LrxMTESLI7hHxT2rXWwBkA/WYFVHUohBCUMSEXR0tLC1paWtKIRQRjDABKHJdWUFAotQ5jDAoK//XIS6tOWTAGpKV9LNeykhD+OpXlNr5lldE+dAZVdnl5+XJ/zNJnq3SybJ9GjVRr7FlyiQnZwcGhzCvj8XgIDg6uUEBCampqACB2aZLwtZqaGnfWW7QOAGRlZXHrUFVVLbZOZmYmtLW1v1qnuLNrQgghRJpKTMiyOOstCx0dHSgqKiI2NlakXPhaT08PfD4fmpqaYl3KycnJyMjI4MaE9fT0xOrk5+cjPj4ePXv25Op8eRnXl9szMTGR2n4RQgghxSkxIR88eLAy4xBTu3ZtmJmZISgoCGPGjOG6kwMDA6GmpgZjY2MAgJWVFS5fvox58+ZxY9iBgYFQVFRE586duTp79+5FamoqGjRoAKDwhiNZWVno2rUrgMJZ1+fOnUNkZCRatiy8NjMiIgKRkZGYNGlSZe46IYSQGqh8g6OVZNKkSbh79y5mzpyJ0NBQbNmyBXv27IGXlxfq1Ckcwxg/fjzevXsHT09PXL58mbspyLBhw7iz/FGjRkFZWRnu7u64ePEijh8/jrlz58LW1hYdO3YEAPTu3Ru6uroYP348zp07h7Nnz2LChAlo3bo1evXqVWVtQAghpGaQeFKXpGPKly5dKncwRVlaWsLHxwfbtm3D5MmToampiXnz5mHcuHFcHX19fezduxfe3t6YNm0aGjZsiLFjx2Lq1KlcHXV1dRw4cACrV6/GnDlzwOfz4ezsjHnz5nF1lJWVsW/fPvz0009YsmQJlJWVYWVlhQULFkBJqcJz3wghhJBSSZxpihtTLigoQFJSEmJiYtCiRQtYWVmVO5DBgwdj8ODBYuWOjo5wdHQsdVkzM7Ov3kdbIBBg//79pdb53//+h+3bt381VkKqg/iVGxEQGlHVYRBCPpM4IZc2pvz48WOMHz+eG7MlhMi/9yPcEPiebp1JiLyQyhiysbExXFxc8PPPP0tjdYQQQkiNI7VJXY0bN0Z0dLS0VkcIkbGGRw+g58PAqg6DEPKZVBLyu3fv4OvrW+XXLhNCJNfsh9mYEvxrVYdBCPmswrOsc3JykJKSgvz8fCxbtkxqgRFCCCE1SYVmWQOFT02ysLBA3759YWdnJ624CCHV0P8a8yv0tKe8vAJkZmZLOSpC5INUZlkTQogkVGorISs7D1EJaWVeVk+7PpSV5PpeRoRUCN3xghBSqaIS0rDo17JfbrV6khUMdBpIPyBC5ITECTk1NRXe3t64fv063r17xz368Es8Hg9Pnz6VaoCEEEJITSBxQl6xYgUuXLiAjh07wsLCAoqKirKMixBCCKlRJE7IN27cgIuLC5YsWSLLeAghleTRy3fl6jomhMiGxDMkatWqxT2WkBBCCCHSJXFCHjRoEE6dOoW8vDxZxkMIIYTUSBJ3WU+fPh1eXl7o2bMnbG1t0ahRI7E6PB4PkydPlmqAhBDZaDXQAZvfZWCmy6aqDoUQgjIk5LNnz+Kvv/5CQUEBfH19i61DCZmQb0edJw/RqqqDIIRwJE7I27dvh46ODhYuXAg9PT2aZU0IIYRIkcQJ+d27d1iwYAG6desmy3gIIYSQGkniSV1t2rRBQkKCLGMhhBBCaiyJE/K8efNw/PhxHDlyBP/++y8KCgpkGRchhBBSo0jcZS18tOLKlSuxcuXKYuvQrTMJIYSQ8pE4IRsZGcHY2FiWsRBCKlHKcFfcfppY1WEQQj6TOCGvXbtWlnEQQipZwqpN+JlunUmI3KCHixJCCCFyoMQz5DZt2sDb2xv9+vUDABgaGoLH45W6MhpDJuTbofL4AfTfRuCVJt0ehBB5UGJCHjhwIHR0dERefy0hE0K+Ha0H9cAWAP1mBVRxJIQQoJSEvGbNGpHXNIZMCCGEyE6FxpDT0tKQnZ0trVgIIYSQGqvUhJybm4ujR49i4cKFIuXh4eHo06cPunTpAlNTU4wfPx6xsbEyDZQQQgipzkpMyDk5ORgzZgyWL1+Os2fPcs9Bjo6OhoeHByIjI2FjYwN3d3dERUVhxIgRSEpKqrTACSGEkOqkxIT8+++/4969e5g7dy7+/vtvKCkVDjf7+PggOzsbffr0wa5duzBv3jycOHECioqK2LFjR6UFTgghhFQnJSbkCxcuoGfPnvDw8ICKigqAwrPmkJAQ8Hg8eHh4cHUbNGiAwYMH48qVKzIPmBBCCKmOSkzIMTExMDMzEym7f/8+Pn78CA0NDbRp00bkPR0dHfz777+yiZIQInUv/YMxY/SGqg6DEPJZiQm5oKAAioqKImV//fUXAKBr165i9dPT01GnTh0ph0cIkZVPxiZ0UxBC5EiJCVlHRwfPnj0TKQsODgaPx4OdnZ1Y/bCwMJEbiRBCCCFEciXeGKRPnz74+eefYWtrCysrK/j5+eHly5do3Lgx7O3tReqePn0a169fx/Tp02UeMCHSxufXhpJS+S7JV1JS/HolOaW9ZBYmP03Ez46TqzoUQghKScju7u64du0apkyZAh6PB8YYatWqhZ9++gnKysoAgIsXL+LQoUO4ffs29PT04O7uXllxEyI1SkoKyMlniEpIK/OybfXUZRBR5VD3OwhngBIyIXKixISsrKyM/fv34/z587h//z74fD769++PVq3+G3N6/Pgx7t69i/79+2PBggXcbGxCvjVRCWlYVI5HEfqu6i2DaAghNVGpz0NWVFREv379uCc+FTVx4kRMnz4dCgr0FEdCCCGkIkpNyF9Ds6oJIYQQ6aBTW0IIIUQOUEImhBBC5ECFuqwJId+uj0btkfAuo6rDIIR8VuIZsq+vL6KjoysxFEJIZYoIuISZLpuqOgxCyGclJmRvb2+Eh4dzrx0cHHDp0qVKCYoQQgipaUq9Djk4OBgdOnRAnTp1kJCQgNevX+P169elrlBLS0vqQRJCCCHVXYkJ+bvvvsOePXsQGhoKAODxeFi9ejVWr15d6gqL3v+aECKf2rXWwBkA/WYFVHUohBCUkpDnzp0Lc3NzvHjxAjk5Ofj555/h6OgIAwODyoyPEEIIqRFKnWVtZ2fHPdnJ398fAwcOhIODQ2XERQghhNQoEl/2FBISAgDIz8/H48ePkZCQAGVlZTRt2hTGxsYyC5AQQgipCcp0HfLly5exYsUKvH37FowxAIVjy02aNMGyZcvEHssoDXl5eejYsSOys7NFyuvWrYt79+4BKHwW8+bNmxEREYFGjRrBxcUF48aNE6n/6NEjeHt74/Hjx+Dz+Rg8eDCmTp2KWrVqcXWio6Oxdu1ahIeHQ1FREc7Ozpg7dy5UVVWlvl+EEELIlyROyOHh4Zg6dSoaNWqEmTNnQl9fH4wxREZG4siRI5g2bRoOHDiAjh07SjXAqKgoZGdnY926dWjRogVXLnygxd27dzFx4kT06tUL06dPx507d+Dt7Q3GGDw8PAAAMTExcHd3h6mpKbZs2YJXr15h8+bNyMjIwNKlSwEAaWlpGDNmDDQ0NLBu3TokJydj/fr1SExMxM6dO6W6T4QQQkhREidkHx8faGtr448//oCamprIe6NGjcKQIUPw66+/Yvfu3VIN8Pnz51BQUEDPnj2LfZjFtm3b0LZtW6xfvx4AYGtri7y8POzYsQOurq5QVlbGrl27oKamhl9++QXKysro1q0bVFRUsGrVKnh5eUFTUxOHDx/Ghw8fEBAQgIYNGwIANDU14enpiQcPHsDExESq+0UIIYR8SeJ7WT98+BBDhw4VS8YAoKqqiu+++w4PHjyQanBA4WVUOjo6xSbj7OxshIeHw8nJSaS8Z8+e+PDhA+7evQsAuH79Orp37w5lZWWujrOzM/Lz8xEWFsbVMTc355IxAFhbW4PP53OXfhFSncSv3IjtPSZVdRgS+19jPpSUFFG/fp1y/+Pza1f1bhBSIqndy5rH4yE3N1daq+O8ePECysrK8PDwwN27d6GkpIRevXph3rx5SExMRG5uLvT09ESW0dXVBVDY3W1iYoI3b96I1VFXV4eqqiqioqIAAJGRkejfv79IHUVFRTRr1oyrU1Y8HlC/vuweUamkpAhAttv4lknaPsJ6Nc37EW4IfH+9qsOQmEptJWRl5yEqIa1cy+tp10ddlVpQUvr6eQh9tkony/bh8aS+ym+GxAnZxMQEf/zxB0aNGoW6deuKvJeRkYHjx4+jXbt2Ug/w+fPnyMjIwNChQzFx4kQ8fvwYPj4+iIqKwqxZswBAbNIVn8/n4kpPTy+2jrBeRkbhzfXT09O/WocQUrWiEtKw6Nfy/YhYPckKRi0bSTkiQqRH4oQ8ZcoUuLm5oW/fvnBxceEmWAkndb19+xYrVqyQeoCbN29G/fr1uRuSmJubo1GjRpg7dy6uXy/8YPJK+EmloKAgMhu8KMYYNzlM0jplwRiQlvaxXMtKQvjrVJbb+JZJ2j419Syo4dED6PkwAoHte1Z1KJUmLy9fos8LfbZKJ8v2adRItcaeJUuckM3MzODj44Mff/wR3t7eXPJijEFDQwObN29Gly5dpB5g586dxcqENysRKnoGK3ytpqbGnfUWd5ablZXFjYmrqqoWWyczMxPa2trlip0Qedbsh9mYAtSohEyIPCvTGLKDgwPs7Ozw5MkTxMfHAwC0tbVhZGQEJSXpP1o5OTkZISEh6NKlC5o3b86Vf/r0CQDQqFEjKCoqIjY2VmQ54Ws9PT3w+XxoamoiJiZGbN0ZGRnc2LKenp5Ynfz8fMTHx6NnT/rCIoQQIltl7otVVFRE+/bt0bt3b/Tu3RsmJiYyScZAYRfy0qVLcejQIZHy8+fPQ1FREV27doWZmRmCgoK4rmkACAwMhJqaGncHMSsrK1y+fBk5OTkidRQVFbkzcCsrK9y6dQupqalcnbCwMGRlZaFr164y2T9CCCFESDaZVErU1dUxevRoHDx4EKqqqjAzM8OdO3ewY8cOjB49Grq6upg0aRLGjh2LmTNnYtCgQbh37x727NmD2bNnc5dKjR8/HufOnYOnpyfGjBmD6OhobNq0CcOGDeMeFzlq1CgcOnQI7u7umDx5MlJTU7F+/XrY2tpK/WYnhBBCSFFynZABYP78+dDU1MSJEyewa9cuaGpqYtq0aRg/fjwAwNLSEj4+Pti2bRsmT54MTU1NzJs3T+TWmfr6+ti7dy+8vb0xbdo0NGzYEGPHjsXUqVO5Ourq6jhw4ABWr16NOXPmgM/nw9nZGfPmzav0fSaEEFLzyH1CrlWrFiZMmIAJEyaUWMfR0RGOjo6lrsfMzAzHjh0rtY5AIMD+/fvLEyYhhBBSIeW7nocQQgghUiVxQnZzc8Nff/3Fvc7IyICbmxuePn0qk8AIIbL16OU79JsVUNVhEEI+K7HL2sbGBkZGRjAyMkLbtm1x+/ZtDBs2jHs/NzcXt2/fRlpa+W5jRwghhJD/lJiQPTw88OzZMwQFBWHnzp3g8Xj48ccfcezYMbRp0wbNmzcHj8cr8S5ZhBBCCJFciQnZ3d2d+/+cnBy0b98ednZ24PP5ePjwIf744w8wxjBx4kS0adMGxsbGaNeundgDGggh8qnVQAdsfpeBmS6bqjoUQggknGUtfGyhjY0N+vXrBwBISUlB165d4eLigvz8fDx58gSnTp2ihEzIN6LOk4doVdVBEEI4JSbkYcOGoU2bNjAyMoKhoSEA0YcvCP/fysoKlpaWMg6TEEIIqd5KTMjm5uZ4/vw5Ll68iJSUFPB4PGzZsgWhoaEwNDSElpYWjSETQgghUlJiQp47dy73/4mJibCzs0Pr1q3x6dMnHD16lHu4xPz582FiYgJjY2MYGxvTfZ8JIYSQcpBoDLlp06YAgN69e3NjyK9fv4a9vT1sbW3x8eNHnDhxAlu2bKHrkgkhhJBykPjWmVpaWqhbty73WlVVFVpaWhg8eDBMTU0BFP/MYUIIIYR8ncQJOSQkROR1vXr1xMpUVVWlExUhROZShrvi9tPEqg6DEPKZ3D9cghBJ8Pm1oaQkeidYJSVFAED9+nVKXVZYr6ZJWLUJP/96varDIIR8RgmZVAtKSgrIyWeISij7rVzb6qnLICJCCCkbSsik2ohKSMOicpzx+a7qLYNo5J/K4wfQfxuBV5p0exBC5AE9fpGQGqr1oB7YcnhOVYdBCPmMEjIhhBAiByghE0IIIXKAEjIhhBAiByghE0IIIXKAEjIhhBAiByghE0IIIXKArkMmpIZ66R+Mn/+4X9VhEEI+o4RMSA31ydgEr67VnAfC/K8xH0pKil+9lSpQ/G1X8/IKkJmZLbP4CKGETAipEVRqKyErO69ct1fV064PZSUa4SOyRQmZkBpKe8ksTH6aiJ8dJ1d1KJWmvLdXXT3JCgY6DaQfECFfoJ98hNRQ6n4H4fzoYlWHQQj5jBIyIYQQIgcoIRNCCCFygBIyIYQQIgcoIRNCCCFygBIyIYQQIgfosidCaqiPRu2R8K7m3BiEEHlHCZmQGioi4FK5rsklhMgGdVkTQgghcoDOkIlc4PNrQ6kCtyYU3nuYEEK+VZSQiVxQUlJATj4r132GAaCtnrqUI6r+2rXWwBkA/WYFVHUohBBQQiZypLz3GQYA31W9pRwNIYRULhpDJoQQQuQAJWRCCCFEDlBCJoQQQuQAjSETQshX/K8xH0pKiqhfv065ls/LK0BmZraUoyLVDSVkQgj5CpXaSsjKzivXVQB62vWhXIFL+kjNQQmZkBoqfuVGBIRGVHUY34zyXgWwepIVDHQaSD8gUu1QQiakhno/wg2B7+nWmYTIC+pHIYQQQuQAJWRCaqiGRw+g58PAqg6DEPIZdVkTqanI/ajpXtSVr9kPszEFQGD7nlUdCiEElJCJFFXkftR0L2pCSE1HCbmIs2fP4tdff0VcXBy0tbXh5eWFgQMHVnVY34zyzkSle1GT6oquYSaSooT8hQsXLmDOnDlwc3ODjY0NgoODMX/+fKioqMDZ2bmqwyOEfIPoGmYiKUrIX9i0aRN69eqFRYsWAQBsbGyQlpaGrVu31piETOPAhEgfXcNMJEEJ+bO4uDjExsZi1qxZIuU9e/bEhQsXEBcXh+bNm1dRdJWHxoEJkR/U3V2z8BhjrKqDkAehoaHw9PTEqVOnYGhoyJU/ffoUgwYNwu7du2Fra1umdX7LTVtQUPbYFRR4VbJsVW77W95nXlrhj658tXqVut2qWLYqt13hvxOPV+blhL7V76CK7PO3jM6QP0tPTwcAqKqqipTz+XwAQEZGRpnX+S0fVIqK5Y+9qpatym1/k/vcoEHh8pW93Spctiq3XdG4y+Nb/g6qiWi2wGfCX5JFD2BhuYICNRUhhBDZoSzzmZqaGgDxM+HMzEyR9wkhhBBZoIT8mZ6eHgAgNjZWpDwmJkbkfUIIIUQWKCF/pquri2bNmuHPP/8UKQ8KCkKLFi2gpaVVRZERQgipCWhS1xcmT56MhQsXon79+rCzs0NISAguXLiAzZs3V3VohBBCqjm67KmIo0ePYu/evXjz5g2aN28OT09PunUmIYQQmaOETAghhMgBGkMmhBBC5AAlZEIIIUQOUEImhBBC5AAlZEIIIUQOUEImhBBC5AAlZDmVmZmJFStWwMrKCqamppgwYQKio6PLtLy9vT1OnTol9l50dDQmTpwIMzMzWFhYYNmyZeV6eEZVKm/7PHr0CK6urjA1NYW1tTU2bdqE3NxckTqLFy+GgYGB2L+iN42RJ2fPnkWfPn3Qvn179OrVCwEBAaXWl6T98vLysGXLFnTr1g0mJiYYNWoUHj58KLudkCFZtE94eHixx4mXl5fsdkRGyto+X1q3bh3c3d3FyqvT8VNpGJFLEyZMYF26dGEnT55kgYGBrF+/fszGxoZ9+PDhq8ump6czV1dXJhAIWEBAgMh7qampzNbWlg0ZMoQFBwczPz8/ZmZmxjw9PWW1KzJRnvaJjo5mHTt2ZB4eHuzKlStsz549zNjYmK1YsUKk3uDBg9n06dPZvXv3RP69f/9exntVPufPn2cGBgbsp59+YlevXmVLly5lAoGAXbhwocRlJGm/5cuXMxMTE3bw4EF26dIl5uLiwkxNTVlsbGxl7JbUyKp9Dh06xDp06CB2nLx69aoydktqytM+QgcPHmQCgYCNGTNG7L3qcvxUJkrIcujvv/9mAoGAhYaGcmXJycmsQ4cObOfOnaUue+3aNebk5MQ6d+5cbEL++eefWYcOHVhKSgpXduXKFSYQCNj9+/eluyMyUt72WbRoEevWrRvLzs7myg4fPszatGnDEhMTGWOM5eXlsfbt27MjR47IbgekrEePHmzGjBkiZdOnT2fOzs7F1pek/eLi4libNm1E2iE7O5vZ2dmxpUuXymAvZEcW7cMYY0uWLGFDhw6VTdCVqKztwxhjiYmJbNasWczQ0JB16tRJLCFXp+OnMlGXtRy6fv06+Hw+rKysuDJ1dXWYm5vj6tWrpS47fvx4GBsbY/fu3SWu29zcHA0bNuTKrK2twefzERoaKp0dkLHyts/169fRvXt3KCsrc2XOzs7Iz89HWFgYACAqKgqfPn2CgYGB7HZAiuLi4hAbGwsnJyeR8p49eyIyMhJxcXFiy0jSfjdv3kR+fj569uzJ1VFWVoadnd1Xj0F5Iqv2AYBnz559M8dJScrTPgCwefNmPH36FPv27UObNm3E3q8ux09lo4QshyIjI6GrqwtFRdFHx+vo6CAqKqrUZU+fPo2NGzdCXV29xHUXfXKVoqIimjVr9tV1y4vytM/Hjx/x5s0bsX1XV1eHqqoqt9zz588BAAEBAbC2toaxsbFcj31FRkYCEH8ama6uLgAU2x6StF9kZCTq168vdhzp6uri9evX+PTpk9T2QZZk1T4FBQV4+fIlEhMTMWjQIBgbG8POzg579+7lnqH+LShP+wCFP/zPnTuHLl26lLje6nD8VDZ6uEQly8vLw7lz50p8v3HjxsjIyICqqqrYe3w+/6uTrwQCQanvp6enl3vdlUFW7ZOeng4AX11OmJDT09OxYcMGfPjwATt37oSbmxuOHTv21fatbCXtF5/PByD+fG9h2dfaobQ6QOGkJxUVlYoFXwlk1T7CnpSoqCjMmjULDRs2xKVLl+Dt7Y2MjAxMmzZN2rsiE+VpHwBo1apVqeutLsdPZaOEXMmys7Mxb968Et/v3LkzatWqVeL7CgoV79Tg8XhiZYwxqay7omTVPsKzlq/t+9ChQ2Fubo5u3bpx73fp0gVOTk7YuXMnNm7cKNF+VJaS9ktYXlx7lHYGJ6xfUp3S2lEeyap9NDU1sXv3brRp0wYaGhoAAEtLS3z69Am7d+/GuHHjik1I8qY87VOW9Uq6PVKIEnIl4/P5ePHiRal1pk2bhvj4eLHyzMzMCn/IVVVVi/3Vm5mZCW1t7QqtWxpk1T7C8uL2PSsrC2pqagAKu9SE3XVC9erVQ8eOHb8aV1UQxl10vzIzM0Xe/5KqqupX209VVZVbR3Hr/RaSDSDb9rG1tRWrY2dnh+PHjyMqKgrt2rWrcPyyVp72kUR1OX4qW9WfEhExenp6iIuLE/uVGRMTIzbWU551x8TEiJTl5+cjPj6+wuuuLOVpHz6fD01NTbF9T05ORkZGBrdcUFBQsZPbsrOzRSbCyQth3LGxsSLlwv0srj0kab+WLVsiNTUVaWlpYnWaNWsmMjFOnsmqfV68eIEjR46IXcMuHBuVx2OlOOVpH0lUl+OnslFClkPW1tb48OEDbty4wZWlpKQgPDwcXbt2rdC6rayscOvWLaSmpnJlYWFhyMrKqvC6K0t528fKygqXL19GTk4OVxYYGAhFRUV07twZAHDy5EksWbJEZNLJ27dvcffuXa6OPNHV1UWzZs3EbloSFBSEFi1aQEtLS2wZSdpP+N/AwECuTk5ODkJDQ7+Z4wSQXfvExMRgxYoVYjOGz58/j2bNmslFb5MkytM+kqgux0+lq9yrrIikXFxcWOfOndmxY8dYUFAQd2OC1NRUrs7Lly/ZkydPil0+Li6u2OuQk5OTmYWFBRswYAALCgpix44dY+bm5mz8+PEy3R9pK0/7REREsHbt2rExY8awkJAQtnfvXmZsbMyWLVvG1bl//z4zMjJi7u7u7MqVK+zMmTPMycmJ2dnZsfT09MrcRYmdOHGCCQQCtmLFChYaGsqWLVvGBAIBO3fuHGOs8G9+7949kfglab/58+ezdu3asX379rGQkBDm6urKTE1NWXR0dKXvY0XIon2ys7PZ4MGDmaWlJTt27Bi7evUqmzt3LjMwMGBBQUFVsp/lVZ72+ZKLi0uxNwapLsdPZaKELKdSU1PZggULmJmZGevYsSObMGGC2B2AXFxcWPfu3YtdvqSEzBhjL168YGPGjGHt27dnlpaW7IcffpDbZFOS8rbP33//zYYOHcqMjY2ZjY0N27hxI8vJyRGr4+Liwjp27MjMzMzYjBkzWEJCgsz3qSJ8fX2Zo6MjMzY2Zr169WL+/v7ce8Iv3Js3b3JlkrRfdnY2++mnn5ilpSUzMTFho0aN+mZuHlOULNonOTmZ/fDDD8zW1pYZGxuzQYMGsYsXL1bWLklVWdvnSyUl5Op0/FQWHmPf0EVzhBBCSDVFY8iEEEKIHKCETAghhMgBSsiEEEKIHKCETAghhMgBSsiEEEKIHKCETAghhMgBSshE6hYsWAADA4Ov/luwYEFVh1puJ0+ehIGBAU6ePFnVoQAAfHx8YGBgUOw9mIWEMd+6dUtq2/3yebm3bt2qtDaJjo7GggULYGtrC2NjY3Tt2hUTJ07E9evXZb7t8qrM9iHfJnq4BJG64cOHw9LSknt9584d+Pn5Yfjw4ejUqRNXrqOjUxXh1Vjm5ubw9vaGvr6+VNbn4eEBDQ0NrF27Virrk9SzZ88watQoqKurY9iwYdDU1ERSUhJOnz6NcePGYcmSJXB1da3UmAiRBkrIROpMTU1hamrKvc7Pz4efnx86dOiAAQMGVGFkNVvz5s3RvHlzqa0vLCwMgwYNktr6JLV+/XrUq1cPAQEBIk8j8vDwwPDhw7Fhwwb07dv3m3nAAyFC1GVNCPmm3Lt3DyYmJmKPBlRWVsbIkSORk5ODZ8+eVVF0hJQfJWRSpXx8fNCuXTtcvHgRVlZWMDU1xfHjx0scEy1a7uPjA1NTU0RERGDs2LHo0KEDbGxssHv3bjDGsGfPHnTv3h0dO3aEh4eHyPoWLFgAR0dH3Lt3D4MHD0b79u3h7OwMX19fqe5jWloaVq5cCRsbGxgbG6NXr174/fffxR7v9+TJE0ydOhVdu3aFkZERLC0tMXv2bCQmJorUi42NxdSpU2Fubg4LCwts3ry5xAfCf6noGLLw9fPnzzF79myYm5vD1NQUkydPLnUsOj4+HgYGBgAAf39/sXHprKwsrFixApaWlujQoQPGjBkj9izpgoIC7N27F87OzjA2NoaNjQ1WrVpV7POqi+Lz+bh16xaioqLE3hsyZAiePHki8kShjIwMbNy4Ec7OzmjXrh1MTU0xbNgwXLp0SWyfTp8+jXXr1qFr164wNTXF999/j5SUFDx8+BAjRoyAiYkJevbsifPnz3PLCseGQ0NDsWDBAnTs2BFdunTBwoULkZKSUuq+VKQdSPVDXdakyuXl5WHJkiXw8PBATk4OOnXqhHPnzkm8fG5uLsaMGYMePXrAyckJJ06cwIYNG3Dz5k0kJCRgzJgxeP/+PX777TcsXLgQBw8e5JZNTU3F+PHj0a1bNwwePBhBQUFYvnw5Pnz4AC8vrwrvW1ZWFlxcXPDmzRuMGjUKTZs2xc2bN7F69WpER0dj2bJlAAqfrztq1Cjo6urC09MTderUwd27d3Hq1Cn8+++/XMxJSUkYMWIEt88qKio4cuQIPnz4UO4YJ02aBH19fcycORNxcXH4/fff8fbtW/zxxx/F1ldXV4e3tzfmzZsHMzMzDBs2DPr6+nj16hUAYMOGDTAwMMDUqVPx9u1b7Nu3Dx4eHggODoaKigoAYPHixQgICMCgQYPg7u6OV69ewdfXF3fv3oWvry9q165dYrxDhgzBjh070LdvX9jY2MDW1hZdunRBy5YtoaioKFKXMQYvLy88ffoULi4u0NHRQWJiIo4ePYqpU6ciMDBQpBt/w4YN0NDQwJQpUxAREYHDhw/j/fv3iIyMxODBg9GvXz8cOHAA8+bNg5GREXR1dbllly9fjrp162LatGl48+YNDh06hMePH+PEiRMlPv+3Iu1AqqGqfLIFqRmET4s5ceKE2Hvbtm1jAoGAbdu2rdjyuLi4UsuFr9euXcvVefnyJRMIBMzU1JQlJydz5bNnz2YGBgYsOzubMVb4eDiBQMBWrVrF1cnLy2OjR49m7du3F3kUYVn2qWi8RkZG7Pnz5yLlGzduZAKBgD179owxxtjSpUuZiYkJe//+vUi9mTNnMoFAwJWvXbuWGRgYsMePH3N1kpKSWJcuXYptr+JiFj61R/h6ypQpIvWWLl3KBAIBi4qKKnXfBAIBmz9/Pvf65s2bTCAQsEGDBrHc3Fyu3MfHhwkEAnbjxg2Rer6+viLru3btGhMIBGz//v2lbjc3N5etXLmSGRoaMoFAwP2zt7dnv/zyC/f3ZazwcZrFbevq1atMIBCwvXv3Msb+ezqara0t+/jxI1dvyJAhTCAQsCNHjnBl169fZwKBgPn5+YnsT7du3USemnbs2DGRbQvrCY+ZirYDqX6oy5rIBWtr6wot36NHD+7/W7RoAQDo2LEj1NXVufJmzZqBMYakpCSRZb88E1ZUVISbmxs+ffok8oD68goKCoJAIICGhgZSUlK4f8J4L1++DKDw7CokJAQNGjTgls3IyODOkLKysgAAV69eRbt27WBkZMTVa9SoEfr06VPuGHv16iXyuk2bNgAg1k6ScnZ2hpLSf51v7dq1E1lfUFAQeDweunXrJtImbdu2hYaGBq5cuVLq+pWUlLBkyRIEBwdj3rx5sLS0hLKyMuLj47FlyxaMHDmSay8TExP8/fffGDx4MLd8fn4+CgoKAACZmZki67axseHO4oH/jiVHR0eurFmzZgCAf//9V2TZUaNGQVVVlXs9aNAg1K9fHyEhIcXuR0XbgVQ/1GVN5EKjRo0qtHzjxo25/xcmg6LrFHZnCr+MAaBBgwYiywLguiETEhIqFBNQON776dMnkcvAvvTmzRsAAI/Hw/v377Fz5068ePECsbGxeP36NTc2LIw5ISEBDg4OYutp2bJluWMsOhtZ2L2an59frvV9+SMIAJfgcnNzARS2CWMMdnZ2xS7P5/Ml2o62tjY8PDzg4eGBT58+4fLly9iyZQseP36MQ4cOwdPTE0Dh8XD06FHcvn0bMTEx3N8EgNjYe9FjprhjSUFBodhlW7VqJbZss2bNSjyOpNUOpPqghEzkgvBL7mtKShJFxw6BwiT3NbVq1RIrEya/4tZZVvn5+ejUqROmTJlS7PtNmjQBAFy5cgXff/89mjRpgi5dunA3vAgLC8POnTu5+jweD9nZ2WLrKZocykLStpfW+goKCsDn87F9+/Zi3y9t3PTOnTsICgrCuHHjoKmpyZWrqKigV69e6NixI+zt7XH37l0AwIcPHzBixAjExcXBysoK9vb2MDQ0hLa2NoYOHSq2/i/P7L9U3mMpPz+/xPaoSDuQ6okSMpFLwi+xnJwckfLydqOWJCkpCZmZmSJnI9HR0QAgMmGnvLS1tZGZmSky6xconHn9119/cdtYuXIldHV1ceLECdStW5erd+bMGZHlmjVrxsX3pS/vmCXvtLW1ERYWBmNjY9SrV0/kvcDAQJFu+6Li4uKwf/9+GBoaFnsNtKamJtTU1Liz8gMHDuDVq1fYv3+/SC+FMGFLU2xsrMjr3NxcJCQklNg7UpF2INUTjSETuaShoQEAeP78OVeWkZGB0NBQqW6HMYbDhw9zr/Py8vD7779DTU2txC/SsrC3t8fz58/FxgN//fVXTJ8+HS9fvgRQONtbS0tLJBm/efMGQUFBAP7rGXBycsLLly9x9epVrl56ejpOnTpV4VjLSkFBQaT7X1L29vYACtvgSyEhIZg2bZrYj5Ciy/L5fPj4+IhdDgYUjsu+f/+e69ZPTU0FINqdzBjDoUOHABT+vaXFz8+P65YHgOPHjyM9PR1OTk7F1q9IO5Dqic6QiVzq0aMHVq1ahR9//BEJCQlQVlbGsWPHRBKWtPzyyy9ISEhA69atceHCBdy7dw8//fQT6tSp89Vl/f39cf/+fbHyNm3aYOTIkfDy8kJQUBCmTJmCESNGoHXr1rhz5w5OnToFW1tb2NraAgBsbW1x/vx5LF26FO3atUN8fDyOHTuGjx8/Avhv8tHYsWNx+vRpTJ06FWPGjIG6ujr8/Pwq1GVdXurq6rh9+zaOHTtWpkl53bp1g4ODA/bu3Yv4+Hh07doVCQkJOHz4MLS0tODh4VHisvXq1cPatWsxa9Ys9O3bF/369YOhoSEKCgoQHh6OCxcuoEePHtwkN1tbWxw8eBBeXl747rvvkJubiwsXLuDx48dQUFAQm9RVEdHR0XBxcUG/fv0QFRUFX19fdO7cucQJdxVpB1I9UUImckldXR27d+/Gxo0bsW3bNjRs2BDDhg1Dy5YtMXPmTKlua8+ePVi+fDn8/f3RqlUrbN++XWRWbWlu376N27dvi5U7ODhg5MiRaNCgAfz8/LBt2zb8+eef8PPzg5aWFr7//nt4enpyXfPCa1hDQkJw6tQpNG3aFAMHDoSjoyNGjhyJmzdvom3btlBVVcWRI0ewfv16+Pn5IT8/H71790br1q2xatUqqbbL18yZMwcbN27EypUrsXLlSvzvf/+TaDkej4etW7fit99+Q0BAAC5fvgx1dXU4OTlh+vTpYpPsinJycsLJkyexb98+XL16FSdPnoSCggJatWqFJUuWYMSIEVy72traYtWqVdi7dy/Wrl2L+vXrw8jICH5+fvjhhx+k+qCNOXPm4P79+9iwYQPU1NTg7u6OadOmlTiGXNF2INUPj1XFT2tC5MCCBQvg7+8vdhcpQsri1q1bcHNzw5o1a0QuryKkrGgMmRBCCJEDlJAJIYQQOUAJmRBCCJEDNIZMCCGEyAE6QyaEEELkACVkQgghRA5QQiaEEELkACVkQgghRA5QQiaEEELkACVkQgghRA78H8ArIyklolPHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(simulations_bias, bins=20)\n",
    "plt.axvline(0, color='red', linestyle='dashed', linewidth=2)\n",
    "plt.title('Approximate Sampling Distribution for a Biased Sample') \n",
    "plt.ylabel('# of Simulations')\n",
    "plt.xlabel('Trump Lead in the Sample')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 341,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.44967"
      ]
     },
     "execution_count": 341,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.count_nonzero(np.array(simulations_bias) > 0) / 100000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the histograms from the two simulations are similar in shape. \n",
    "They are symmetric with reasonable length tails, i.e., they appear to roughly follow the normal curve.\n",
    "The second histogram is shifted slightly to the left, which reflects the non-response bias we introduced."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "q2c",
     "locked": true,
     "schema_version": 2,
     "solution": false
    }
   },
   "source": [
    "##  Would increasing the sample size have helped?\n",
    "\n",
    "With our simulation study we can get insight into the answer to this question. \n",
    "For example, we can try a sample size of 12,000 and run 100,000 simulations for both scenarios."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 342,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulations_big = [trump_advantage(votes, 12000) for i in range(100000)] \n",
    "simulations_bias_big = [trump_advantage(votes_bias, 12000) for i in range(100000)] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 343,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.78798 0.3711\n"
     ]
    }
   ],
   "source": [
    "scenario_no_bias = np.count_nonzero(np.array(simulations_big) > 0) / 100000\n",
    "scenario_bias = np.count_nonzero(np.array(simulations_bias_big) > 0) / 100000\n",
    "print(scenario_no_bias, scenario_bias)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By analyzing over 4,000 polls for 600 state-level, gubernatorial, senatorial, and presidential elections, {cite}`shirani2018` et al found that on average these polls exhibited a bias of about 1.5 percentage points.\n",
    "\n",
    "When the margin of victory is relatively small as it was in 2016, a larger sample size reduces the sampling error, but unfortunately, if there is bias, then the predictions are close to the biased estimate. If the bias pushes the prediction from one candidate (Trump) to another (Clinton), then we have a \"surprise\" upset. "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
